// ------------------------------------------------------------------- //
//              Canadian Catch Age Model (ccam)                        //
//                        a port of                                   //
//         integrated Statistical Catch Age Model (iscam)              //
//                                                                     //
//                           VERSION 1.1                               //
//               Tue Jul 19 22:23:58 PDT 2011                          //
//                                                                     //
//                                                                     //
//           Created by Steven Martell on 2010-04-09                   //
//           Copyright (c) 2010. All rights reserved.                  //
//                                                                     //
// AUTHORS: SJDM Steven Martell                                        //
//		      RF Robyn Forrest                                           //
//		      CG Chris Grandin                                           //
//                                                                     //
// CONVENTIONS: Formatting conventions are based on the The            //
//               Elements of C++ Style (Misfeldt et al. 2004)          //
//                                                                     //
// NAMING CONVENTIONS:                                                 //
//             Macros       -> UPPERCASE                               //
//             Constants    -> UpperCamelCase                          //
//             Functions    -> lowerCamelCase                          //
//             Variables    -> lowercase                               //
//                                                                     //
// CHANGED   add option for using empirical weight-at-age data         //
// TODO:     add gtg options for length based fisheries                //
// CHANGED   add time varying natural mortality rate with splines      //
// TODO:     add cubic spline interpolation for tm varying M           //
// CHANGED   fix the type 6 selectivity implementation. not working.   //
// TODO:     fix cubic spline selectivity for only years when data avai//
// CHANGED:  fixed a bug in the simulation model log_ft_pars goes out  //
//           of bounds.                                                //
// TODO:     write a projection routine and verify equilibrium calcs   //
// TODO:     add DIC calculation for MCMC routines (in -mcveal phase)  //
//                                                                     //
//                                                                     //
//                                                                     //
//                                                                     //
//                                                                     //
// ------------------------------------------------------------------- //
//-- CHANGE LOG:                                                     --//
//--  Nov 30, 2010 -modified survey biomass by the fraction of total --//
//--                mortality that occurred during the time of the   --//
//--                survey. User specifies this fraction (0-1) in the--//
//--                data file as the last column of the relative     --//
//--                abundance index.                                 --//
//--                                                                 --//
//--  Dec 6, 2010 -modified the code to allow for empiracle weight-  --//
//--               at-age data to be used.                           --//
//--              -rescaled catch and relative abundance /1000, this --//
//--               should be done in the data file and not here.     --//
//--                                                                 --//
//--  Dec 20, 2010-added prior to survey q's in control file         --//
//--                                                                 --//
//--  Dec 24, 2010-added random walk for natural mortality.          --//
//--                                                                 --//
//--  Jan 23, 2011-in Penticton Hospital with my mom in ICU, adopting--//
//--               the naming conventions in The Elements of C++     --//
//--               style to keep my mind busy.                       --//
//--                                                                 --//
//-- May 5, 2011- added logistic selectcitivty as a fucntion of      --//
//--              mean body weight.  3 parameter logistic.           --//
//--              NOT WORKING YET                                    --//
//--                                                                 --//
//-- May 6, 2011- added pre-processor commands to determin PLATFORM  --//
//--              either "Windows" or "Linux"                        --//
//--                                                                 --//
//-- use -mcmult 1.5 for MCMC with log_m_nodes with SOG herrning     --//
//--                                                                 --//
//-- July 2012 - used for pacific cod exploratory assessment work - after lots of modifications Oct 2011 - March 2012 for pacific hake assessment (annotated inline)
// ------------------------------------------------------------------- //


DATA_SECTION
	init_adstring DataFile;
	init_adstring ControlFile;

  !! BaseFileName = adstring("ccam");
	!! ReportFileName = BaseFileName + adstring(".rep");
  !! ad_comm::change_datafile_name(DataFile);

	int SimFlag;
	int rseed;
	int retro_yrs;
	int useManagementParameters;  // use the management-based model or
																// population-based model.  Set this
																// value in the control file (misc variable 14)

	LOC_CALCS
		SimFlag=0;
		rseed=999;
		int on,opt;
		//the following line checks for the "-SimFlag" command line option
		//if it exists the if statement retreives the random number seed
		//that is required for the simulation model
		if((on=option_match(ad_comm::argc,ad_comm::argv,"-sim",opt))>-1)
		{
      cout<<"-sim flag detected. argv[2] = "<<argv[2]<<"\n";
			SimFlag=1;
			rseed=atoi(ad_comm::argv[on+1]);
      cout<<"Simulation seeded.\n";
			//if(SimFlag)exit(1);
		}
		// command line option for retrospective analysis. "-retro retro_yrs"
		retro_yrs=0;
		if((on=option_match(ad_comm::argc,ad_comm::argv,"-retro",opt))>-1)
		{
			retro_yrs=atoi(ad_comm::argv[on+1]);
			cout<<"______________________________________________________\n"<<endl;
			cout<<"    **Implementing Retrospective analysis** "<<endl;
			cout<<"    **Number of retrospective years = "<<retro_yrs<<endl;
			cout<<"______________________________________________________"<<endl;
		}
	END_CALCS

	//Read in objects from data file using the init_ prefix
	init_int syr;
	init_int nyr;
	!! cout<<"syr\t"<<syr<<endl;
	!! cout<<"nyr\t"<<nyr<<endl;

	init_int sage;
	init_int nage;
	vector age(sage,nage);
	!! age.fill_seqadd(sage,1);

	init_int ngear;				//number of gear types with unique selectivities
	!! cout<<"ngear\t"<<ngear<<endl;
	init_vector allocation(1,ngear);
	ivector fsh_flag(1,ngear);
	LOC_CALCS
		//If allocation >0 then set fish flag =1 else 0
		int k;
		allocation = allocation/sum(allocation);
		for(k=1;k<=ngear;k++)
		{
			if(allocation(k)>0)
				fsh_flag(k)=1;
			else
				fsh_flag(k)=0;
		}
	END_CALCS

	init_number fixed_m;
	init_number linf;
	init_number vonbk;
	init_number to;
	init_number a;
	init_number b;
	init_number ah;
	init_number gh;

	vector la(sage,nage);		//length-at-age
	vector wa(sage,nage);		//weight-at-age

	LOC_CALCS
	  cout<<"linf\t"<<linf<<endl;
	  la = linf*(1.-exp(-vonbk*(age-to)));
	  wa = a*pow(la,b);
	  cout<<setprecision(2);		//2 decimal places for output
	  cout<<"la\n"<<la<<endl;
	  cout<<"wa\n"<<wa<<endl;
	  cout<<setprecision(7);
	END_CALCS

	//Time series data
	init_matrix catch_data(syr,nyr,1,ngear+1);
	matrix obs_ct(1,ngear,syr,nyr);
	int ft_count;
	int i;
	LOC_CALCS
		ft_count=0;
		for(k=1;k<=ngear;k++)
			obs_ct(k)=column(catch_data,k+1);
		for(k=1;k<=ngear;k++)
		{
			for(i=syr;i<=nyr;i++)
				if(obs_ct(k,i)>0) ft_count++;
		}
		cout<<"ft_count\n"<<ft_count<<endl;
		cout<<"last row of catch \n"<<catch_data(nyr)<<endl;
		cout<<"Ok after catch extraction"<<endl;
	END_CALCS

	init_int nit;
	init_ivector nit_nobs(1,nit);
	//#survey type
	//## 1 = survey is proportional to vulnerable numbers
	//## 2 = survey is proportional to vulnerable biomass
	//## 3 = survey is proportional to spawning biomass (e.g., herring spawn survey)
	init_ivector survey_type(1,nit);
	init_3darray survey_data(1,nit,1,nit_nobs,1,5);
	//init_matrix survey_data(1,nit,1,4);
	imatrix iyr(1,nit,1,nit_nobs);
	imatrix igr(1,nit,1,nit_nobs);
	matrix it(1,nit,1,nit_nobs);
	matrix it_wt(1,nit,1,nit_nobs);		//relative weight
	matrix it_timing(1,nit,1,nit_nobs);	//timing of the survey (0-1)
	LOC_CALCS
		for(i=1;i<=nit;i++)
		{
			iyr(i)=ivector(column(survey_data(i),1));
			igr(i)=ivector(column(survey_data(i),3));
			it(i)=column(survey_data(i),2);
			it_wt(i)=column(survey_data(i),4)+1.e-10;//add a small constant to allow 0 weight
			it_timing(i)=column(survey_data(i),5);
		}
		cout<<"Last row of the relative abundance data\n"<<survey_data(nit)(nit_nobs(nit))<<endl;
		cout<<"OK after relative abundance index"<<endl;

    cout<<"Survey data:"<<endl;
    cout<<survey_data<<endl;
    cout<<"survey_type:"<<endl<<survey_type<<endl;
    cout<<"igr:"<<endl<<igr<<endl;

	END_CALCS

	// Age comps for all gears.
	init_int na_gears;	//total number of aging observations
	init_ivector na_nobs(1,na_gears);
	init_ivector a_sage(1,na_gears);	//youngest age in the ageing matrix
	init_ivector a_nage(1,na_gears);	//oldest age in the ageing matrix
	
	init_3darray A(1,na_gears,1,na_nobs,a_sage-2,a_nage);

	//Mean weight-at-age data (units are kg) (if exists)
	init_int n_wt_nobs; //read in from data file
	init_matrix tmp_wt_obs(1,n_wt_nobs,sage-1,nage); //read in from data file
  //!! cout<<"**********************Creating wt_obs, nyr = "<<nyr<<" retro_yrs = "<<retro_yrs<<endl;
	matrix wt_obs(syr,nyr+1-retro_yrs,sage,nage);		//weight-at-age
	matrix wt_dev(syr,nyr+1-retro_yrs,sage,nage);		//standardized deviations in weight-at-age
	matrix fec(syr,nyr+1-retro_yrs,sage,nage);		       //fecundity-at-age
	vector avg_fec(sage,nage);				              //average fecundity-at-age
	vector avg_wt(sage,nage);				                // average weight at age (for management oriented mode) 
                                                  // - to be consistent with use of average fecundity

	LOC_CALCS
		int iyr;
		avg_fec.initialize();
		//if no weight at age data
		for(i=syr;i<=nyr+1-retro_yrs;i++)
		{
			wt_obs(i)=wa;			//wa is weight at age from length-weight relationship created above
			fec(i)=elem_prod(plogis(age,ah,gh),wt_obs(i));
		}

		//if empirical weight-at-age data exist, then overwrite wt_obs & fec.
    for(i=1;i<=(n_wt_nobs-retro_yrs);i++)
		{
			iyr=tmp_wt_obs(i,sage-1);  //index for year - get from data file
			wt_obs(iyr)=tmp_wt_obs(i)(sage,nage); // overwrite wt_obs with data from data file
			fec(iyr)=elem_prod(plogis(age,ah,gh),wt_obs(iyr));
			cout<<"i = "<<i<<endl;
		}

		int nfec = fec.rowmax()-fec.rowmin()+1;
		int nw = wt_obs.rowmax()-wt_obs.rowmin()+1;

		//from Jake Schweigert: use mean-weight-at-age data
		//from the last 5 years for the projected mean wt.
		dvector tmp=colsum(wt_obs.sub(nyr-5-retro_yrs,nyr-retro_yrs))/6.;
		wt_obs(nyr+1-retro_yrs) = tmp;
    fec(nyr-retro_yrs+1)=elem_prod(plogis(age,ah,gh),wt_obs(nyr-retro_yrs+1)); //set fecundity in 2011

    //GET AVERAGE WEIGHT OVER TIME SERIES FOR EQUILIBRIUM CALCULATIONS (do not include nyr+1)
		avg_wt=colsum(wt_obs.sub(syr,nyr-retro_yrs))/nw;
		avg_fec=colsum(fec.sub(syr,nyr-retro_yrs))/(nfec-1);
		
		//July 14, 2011 Error handler to ensure non-zero wt_obs
		//in the data. Otherwise causes and error with weight-based
		//selectivities.
		if(min(wt_obs)==0)
		{
			cout<<"Cannot have a observed 0 mean weight at age\n";
			cout<<"in the data file.  Please fix.\n Aborting program!"<<endl;
			exit(2);
		}

		//May 5, 2011 SJDM: Calculating standardized deviates
		//in mean weights-at-age from empiracle data to use
		//with selectivity as a function of mean weight at age.
		//Idea borrowed from Vivian Haist in HCAM model.
		/*
			CHANGED: selectivity function based on average weight.
			Based on the logistic function;
			1) calculate a matrix of standardized deviates
			wt_dev = (wt_obs-mean(wt_obs))/sd(wt_obs);

			Not implemented, using the version provided by Vivian.

			Aug 5, 2011, noticed that Vivian's method was implemented
			in an awkward way where GN selectivities were a random walk
			sel(t+1)=sel(t) + tmp(2);

			Trying a compromize where estimating an additional parameter
			that attemps to explain residual variation in age-comps
			via changes in standardized mean weights at age. This is 
			implemented as:

			log_sel = log(plogis(age,s1,s2)) + s3*delta
			where delta is a matrix of standardized deviations
			(mu=0, std=1) of weights at age.  delta is calculated
			based on the wt_dev matrix above.
		*/
		wt_dev.initialize();
		dmatrix mtmp = trans(wt_obs);
    cout<<endl<<mtmp<<endl;

		for(i=sage;i<=nage;i++)
		{
			dvector wa_dev = (mtmp(i)-mean(mtmp(i)))/sqrt(var(mtmp(i)));
			mtmp(i) = wa_dev;
		}
		wt_dev = trans(mtmp);	//each column has mean=0 sd=1

	END_CALCS

	init_int eof;
	LOC_CALCS
	  cout<<"eof = "<<eof<<endl;
	  if(eof==999){
		cout<<"\n -- END OF DATA SECTION -- \n"<<endl;
	  }else{
		cout<<"\n *** ERROR READING DATA *** \n"<<endl; exit(1);
	  }
	END_CALCS

	// fmsy and msy moved to parameter section - use value(fmsy) in call to equilibrium in calc_reference_points
	//number fmsy;					// Fishing mortality rate at Fmsy
	//number msy;						// Maximum sustainable yield
	// End change

	number bmsy;					//Spawning biomass at MSY
	vector age_tau2(1,na_gears);	//MLE estimate of the variance for the age comps
	3darray d3C(1,ngear,syr,nyr,sage,nage);		//catch-age for simulation model (could be declared locally 3d_array)

	// ***************************************************
	// ** Read parameter controls from ControlFile
	// ***************************************************
	!! ad_comm::change_datafile_name(ControlFile);

	init_int npar;
	init_matrix theta_control(1,npar,1,7);
	!! cout<<theta_control<<endl;

	vector theta_ival(1,npar);
	vector theta_lb(1,npar);
	vector theta_ub(1,npar);
	ivector theta_phz(1,npar);
	ivector theta_prior(1,npar);
	LOC_CALCS
		theta_ival = column(theta_control,1);
		theta_lb = column(theta_control,2);
		theta_ub = column(theta_control,3);
		theta_phz = ivector(column(theta_control,4));
		theta_prior=ivector(column(theta_control,5));
	END_CALCS

	// ***************************************************
	// ** Read selectivity parameter options
	// ***************************************************
	// type 1 = logistic (2pars)
	// type 2 = selcoffs (A-1 pars)
	// type 3 = cubic spline (age_nodes)
	// type 4 = time varying cubic spline (age_nodes for each year)
	// type 5 = bicubic spline with age_nodes adn yr_nodes
	// type 6 = fixed logistic by turning sel_phz to (-ve)
	// type 7 = logistic (3pars) as a function of body weight.
	init_ivector isel_type(1,ngear);  	//Switch for selectivity
	ivector isel_npar(1,ngear);			//ivector for # of parameters for each gear.
	ivector jsel_npar(1,ngear);			//ivector for the number of rows for time-varying selectivity.
	init_vector ahat(1,ngear);			//age-at-50% vulnerbality for logistic function
	init_vector ghat(1,ngear);			//std at 50% age of vulnerability for logistic funciton
	init_vector age_nodes(1,ngear);		//No. of age-nodes for bicubic spline.
	init_vector yr_nodes(1,ngear);		//No. of year-nodes for bicubic spline.
	init_ivector sel_phz(1,ngear);		//Phase for estimating selectivity parameters.
	init_vector sel_2nd_diff_wt(1,ngear);	//Penalty weight for 2nd difference in selectivity.
	init_vector sel_dome_wt(1,ngear);		//Penalty weight for dome-shaped selectivity.
	!! cout<<"isel_type\n"<<isel_type<<endl;
	!! cout<<"ahat\n"<<ahat<<endl;
	!! cout<<"ghat\n"<<ghat<<endl;
	!! cout<<"age_nodes\n"<<age_nodes<<endl;
	!! cout<<"yr_nodes\n"<<yr_nodes<<endl;
	!! cout<<"est phase\n"<<sel_phz<<endl;
	!! cout<<"sel_2nd_diff_wt\n"<<sel_2nd_diff_wt<<endl;
	!! cout<<"sel_dome_wt\n"<<sel_dome_wt<<endl;
	
	LOC_CALCS
		//cout up the number of selectivity parameters
		//depending on the value of isel_type
		isel_npar.initialize();
		for(i=1;i<=ngear;i++)
		{
			jsel_npar(i)=1;
			switch(isel_type(i))
			{
				case 1:
					// logistic selectivity
					isel_npar(i) = 2;
					break;
				case 2:
					// age-specific coefficients
					isel_npar(i) = (nage-sage);
					break;
				case 3:
				 	// cubic spline
					isel_npar(i) = age_nodes(i);
					break;
				case 4:
					// annual cubic splines
					isel_npar(i) = age_nodes(i);
					jsel_npar(i) = (nyr-syr-retro_yrs)+1;
					break;
				case 5:
					// bicubic spline
					jsel_npar(i) = age_nodes(i);
					isel_npar(i) = yr_nodes(i);
					break;
				case 6:
					// fixed logistic (no parameters estimated)
					// ensure sel_phz is set to negative value.
					isel_npar(i) = 2;
					if(sel_phz(i)>0) sel_phz(i) = -1;
					break;
				case 7:
					// CHANGED: Now working, Vivian Haist fixed it.
					// logistic (3 parameters) with mean body
					// weight deviations.
					isel_npar(i) = 2;
					break;
				case 8:
					// Alternative logistic selectivity with wt_dev coefficients.
					isel_npar(i) = 3;
					break;
				default: break;
			}
		}
				
	END_CALCS

	//RF CG change
	//control parameters for GAMMA prior on SURVEY ghat (std of logistic)
	init_number ghat_p1; //parameter 1 of gamma prior
	init_number ghat_p2; //parameter 2 of gamma prior
	init_vector ghat_pswitch(1,ngear); //switch 1=use prior; 2=don't use prior - for each gear
	!! cout<<"ghat_p1\n"<<ghat_p1<<endl;
	!! cout<<"ghat_p2\n"<<ghat_p2<<endl;
	!! cout<<"ghat_pswitch\n"<<ghat_pswitch<<endl;
	// Controls for prior on survey q.
	init_int nits;
	init_ivector q_prior(1,nits);
	init_vector q_mu(1,nits);
	init_vector q_sd(1,nits);
	!! cout<<"nits\n"<<nits<<endl;
	!! cout<<"q Prior\n"<<q_mu<<endl<<q_sd<<endl;

	// Miscellaneous controls (cntrl vector)
	// 1 -> verbose
	// 2 -> recruitment model (1=beverton-holt, 2=rickers)
	// 3 -> std in catch first phase
	// 4 -> std in catch in last phase
	// 5 -> assumed unfished in first year (0=FALSE, 1=TRUE)
	// 6 -> minimum proportion at age to consider in the dmvlogistic likelihood
	// 7 -> mean fishing mortality rate to regularize the solution
	// 8 -> standard deviation of mean F penalty in first phases
	// 9 -> standard deviation of mean F penalty in last phase.
	// 10-> phase for estimating deviations in natural mortality.
	// 11-> std in natural mortality deviations.
	// 12 number of estimated nodes for deviations in natural mortality	
	// 13-> fraction of total mortality that takes place prior to spawning
	// 14-> switch for age-composition likelihood (1=dmvlogistic,2=dmultinom)
	// 15-> Use management parameters (1) or population parameters (0).

	init_vector cntrl(1,15);
	init_int pyrs; //number of projection years
	init_int hcr; //harvest control rule index 1. 40-10 Rule (only rule implemented so far)
	int pyr;

	!!pyr = nyr+pyrs;
	init_vector SSstream(nyr+1,pyr); //catch stream from SS
	init_int eofc;
	int verbose;

	LOC_CALCS
		verbose = cntrl(1);

		useManagementParameters = cntrl(15);
		if(useManagementParameters) {
			cout<<"******* MANAGEMENT ORIENTED MODE *******"<<endl;
			cout<<"                                                                           "<<endl;
			if (cntrl(2) != 1) cout<<"********** WARNING: Only Beverton-Holt available for management-oriented case. Ctrl(2) has been changed to 1. ************"<<endl;
			cntrl(2)=1;
		}

		cout<<"cntrl\n"<<cntrl<<endl;
		cout<<"eofc\t"<<eofc<<endl;
		if(eofc==999){
			cout<<"\n -- END OF CONTROL FILE -- \n"<<endl;
		}else{
			cout<<"\n ***** ERROR CONTROL FILE ***** \n"<<endl; exit(1);
		}
	END_CALCS

	int nf;
	ivector ilvec(1,6);
	!!ilvec=ngear;			//number of fisheries
	!!ilvec(2)=nit;			//number of surveys
	!!ilvec(3)=na_gears;	//number of age-comps
	!!ilvec(4)=1;
	// SM Oct 31, 2010.  Implementing retrospective analysis.
	!! nyr = nyr - retro_yrs;

	//------------REFERENCE POINTS------------------------------
	number sb40_sb;		// spawning biomass at sb40 reference level
	number sb40_spr;	// spr value at sb40 reference level
	number sb40_u;		// exploitation fraction (catch/3+biomass)
	number sb40_y;		// yield at sb40 reference level
	number sb40_de;		// Test - spawning depletion at sb40 reference level

	number spr40_sb;
	number spr40_spr;
	number spr40_u;
	number spr40_y;

	number msy_sb;		// spawning biomass at msy
	number msy_spr;		// spr at msy
	number msy_u;		  // exploitation fraction at msy
	number msy_y;		  // maximum sustainable yield.

PARAMETER_SECTION
	//Leading parameters
	//theta[1]		log_ro, or log_msy
	//theta[2]		steepness(h), or log_fmsy
	//theta[3]		log_m, or m
	//theta[4]		log_avgrec
  //theta[5]		log_recinit
	//theta[6]		rho
	//theta[7]		varphi (precision - was called kappa but kappa is recruitment compensation ratio in this tpl)

	init_bounded_number_vector theta(1,npar,theta_lb,theta_ub,theta_phz);
	!! for(int i=1;i<=npar;i++) theta(i)=theta_ival(i);
	//Selectivity parameters (A very complicated ragged array)
	init_bounded_matrix_vector sel_par(1,ngear,1,jsel_npar,1,isel_npar,-15.,15.,sel_phz);

	LOC_CALCS
		//initial values for logistic selectivity parameters
		//set phase to -1 for fixed selectivity.
		for(int k=1;k<=ngear;k++)
		{
			if( isel_type(k)==1 || isel_type(k)==6 || isel_type(k)>=7 )
			{
				sel_par(k,1,1) = log(ahat(k));
				sel_par(k,1,2) = log(ghat(k));
			}
			//cout<<sel_par(k)<<endl;
		}
		//exit(1);
	END_CALCS

	//Fishing mortality rate parameters
	init_bounded_vector log_ft_pars(1,ft_count,-30.,3.0,1);

	LOC_CALCS
		if(!SimFlag) log_ft_pars = log(0.1);
	END_CALCS

	//CHANGED Vivian Haist: possible bias in Bo due to dev_vector for log_rec_devs
	//	-Try and compare estimates using init_bounded_vector version.
	//	-Now estimating Rinit, Rbar and Ro.
	// May 16, Just looked at Ianelli's code, he used bounded_vector, and
	// separates the problem into init_log_rec_devs and log_rec_devs

	//Annual recruitment deviations
	//!! int ii;
	//!! ii=syr-nage+sage;
	//!! if(cntrl(5)) ii = syr;  //if initializing with ro
	//init_bounded_vector log_rec_devs(ii,nyr,-15.,15.,2);

	!! int init_dev_phz = 2;
	!! if(cntrl(5)) init_dev_phz = -1;
	init_bounded_vector init_log_rec_devs(sage+1,nage,-15.,15.,init_dev_phz);
	init_bounded_vector log_rec_devs(syr,nyr,-15.,15.,2);

	//Deviations for natural mortality
	!! int m_dev_phz = -1;
	!! m_dev_phz = cntrl(10);
	!! int n_m_devs = cntrl(12);
	init_bounded_vector log_m_nodes(1,n_m_devs,-5.0,5.0,m_dev_phz);
	//init_bounded_vector log_m_devs(syr+1,nyr,-5.0,5.0,m_dev_phz);

	objective_function_value f;

  number ro;					//unfished age-1 recruits
	number bo;					//unfished spawning stock biomass
	number kappa;
	number m;					//initial natural mortality rate
	number m_bar;				//average natural mortality rate
	number log_avgrec;			//log of average recruitment.
	number log_recinit;			//log of initial recruitment in syr.
	//number log_avg_f;			//log of average fishing mortality DEPRICATED
	number rho;					//proportion of the observation error
	number varphi				//total precision in the CPUE & Rec anomalies.
	number so;
	number beta;
	vector log_rt(syr-nage+sage,nyr);
	vector vax(sage,nage);		//survey selectivity coefficients
	vector q(1,nit);			//survey catchability coefficients
	vector sbt(syr,nyr+1);		//spawning stock biomass
	vector rt(syr+sage,nyr); 	//predicted sage recruits from S-R curve
	vector delta(syr+sage,nyr);	//residuals for stock recruitment
	vector avg_log_sel(1,ngear);//conditional penalty for objective function
	vector log_m_devs(syr+1,nyr);// log deviations in natural mortality
	matrix nlvec(1,6,1,ilvec);	//ragged matrix for negative loglikelihoods
	//matrix jlog_sel(1,ngear,sage,nage);		//selectivity coefficients for each gear type.
	//matrix log_sur_sel(syr,nyr,sage,nage);	//selectivity coefficients for survey.
	matrix N(syr,nyr+1,sage,nage);			//Numbers at age
	matrix F(syr,nyr,sage,nage);			//Age-specific fishing mortality
	matrix M_tot(syr,nyr,sage,nage);		//Age-specific natural mortality
	matrix ft(1,ngear,syr,nyr);				//Gear specific fishing mortality rates
	matrix log_ft(1,ngear,syr,nyr);			//Gear specific log fishing mortlity rates
	matrix Z(syr,nyr,sage,nage);
	matrix S(syr,nyr,sage,nage);
	matrix ct(1,ngear,syr,nyr);				//predicted catch biomass
	matrix epsilon(1,nit,1,nit_nobs);		//residuals for survey abundance index
	matrix pit(1,nit,1,nit_nobs);			//predicted relative abundance index
	matrix qt(1,nit,1,nit_nobs);			//catchability coefficients (time-varying)
	3darray Ahat(1,na_gears,1,na_nobs,a_sage-2,a_nage);		//predicted age proportions by gear & year
	3darray A_nu(1,na_gears,1,na_nobs,a_sage-2,a_nage);		//residuals for age proportions by gear & year
	3darray log_sel(1,ngear,syr,nyr,sage,nage);		//selectivity coefficients for each gear type.
	3darray Chat(1,ngear,syr,nyr,sage,nage);		//predicted catch-at-age
	//Miscellaneous moved or extra or moved parameters
	// fmsy and msy moved to parameter section - use value(fmsy) in call to equilibrium in calc_reference_points
	number fmsy;					//Fishing mortality rate at Fmsy
	number msy;					//Maximum sustainable yield
	number h; //not sure if this is needed or available elsewhere (due to if statement)
	number Pf; //priors component of objective function
	number Lf;  //likelihood and penalties component of objective function
	number f40; //Fishing mortality that results in SPR=0.4
	number fb40; //Fishing mortality that results in SB/SB0=0.4
	number spr_msy; //spr at msy
	vector ispr(syr,nyr); //SPR under history of fishery (ignore rows that relate to survey)
	vector priors(1,npar+1); //MOVED FROM CALC_OBJECTIVE FUNCTION +1 to allow for survey ghat prior
	vector qvec(1,nits); //MOVED FROM CALC_OBJECTIVE FUNCTION
	vector pvec(1,7); //MOVED FROM CALC_OBJECTIVE FUNCTION
	sdreport_number sd_depletion;

PRELIMINARY_CALCS_SECTION
  //Run the model with input parameters to simulate real data.
  nf=0;
  if(SimFlag)
  {
     calcSelectivities(); //need to calculate selectivity first
     if(useManagementParameters)
     {
        initParameters_M();
     }else{
        initParameters();
     }
     calcTotalMortality();
     simulation_model(rseed);
  }

RUNTIME_SECTION
    maximum_function_evaluations 100,200,500,25000,25000
    convergence_criteria 0.01,0.01,1.e-5,1.e-5

PROCEDURE_SECTION
	//CG RF change - calc selectivities called before init model (selectivity needed for management oriented approach)
	calcSelectivities();
	
	if(useManagementParameters){
			initParameters_M();
		}else{
			initParameters();
		}
	
	calcTotalMortality();
	calcNumbersAtAge();
	calcFisheryObservations();
	calcAgeProportions();
	
	calc_survey_observations();
	calc_stock_recruitment();
	
	calc_objective_function();
	sd_depletion=sbt(nyr)/bo;

	if(mc_phase())
	{
		mcmcPhase=1;
	}

	if(mceval_phase())
	{
		mcmcEvalPhase=1;
		mcmc_output();
	}

	//The following causes a linker error
	//duplicate symbol in libdf1b2o.a
	//dvariable a=3.0;
	//cout<<"testing gammln(dvariable)"<<gammln(a)<<endl;

  //initParameters for management-oriented case
  //derives ro and steepness from leading management parameters
  //Calls function init_model_BH in managementoriented.cxx
FUNCTION initParameters_M
	msy=theta(1); //msy, fmsy and m estimated in normal space with lognormal prior
	fmsy=theta(2);
	m = mfexp(theta(3));
	// m = theta(3); //RF CG change for testing against TINSS - estimate M directly
	log_avgrec = theta(4);
	log_recinit = theta(5);
	rho=theta(6);
	varphi=theta(7); //RF NOTE: THIS WAS INCORRECTLY CALLED KAPPA IN CTRL FILE

 	int fisheryFlag;
 	for(fisheryFlag=1;fisheryFlag<=ngear;fisheryFlag++){
 		if(fsh_flag(fisheryFlag))
      break;
 	}
 	// NOTE - LOOP ABOVE GETS INDEX "fisheryFlag" TO SET LOG_SELECTIVITY FROM FIRST COMMERCIAL GEAR TYPE
 	// Need to add catch from other commercial fisheries together and add warning flag

 	// CURRENTLY USES SELECTIVITY FROM NYR, AND AVERAGE WEIGHT AND FECUNDITY
 	// if option 1 or 6, selectivity is fixed through time anyway
 	init_model_BH(msy,fmsy,ro,kappa,m,age,avg_wt,avg_fec,log_sel(fisheryFlag)(nyr)); //function returns ro and kappa as function of fmsy and msy - BEVERTON HOLT CASE ONLY
  //init_model_BH(msy,fmsy,ro,kappa,m,age,wt_obs(syr),fec(syr),log_sel(fisheryFlag)(nyr)); //function returns ro and kappa as function of fmsy and msy - BEVERTON HOLT CASE ONLY

  h = kappa/(kappa + 4.);

FUNCTION initParameters
  {
	  /*
	  This function is used to extract the specific parameter values
	  from the init_bounded_number_vector to the specific variables
	  used in the code.

	  Note that you must call this routine before runnning the
	  simulation model to generate fake data.
	  */

	  ro = mfexp(theta(1));
	  h = theta(2);
	  switch(int(cntrl(2)))
	  {
		  case 1:
			  //Beverton-Holt model
			  kappa = (4.*h/(1.-h));
			  break;
		  case 2:
			  //Ricker model
			  kappa = pow((5.*h),1.25);
		  break;
	  }

	  // TODO Alternative parameterization using MSY and FMSY as leading parameters
	  m = mfexp(theta(3));
	  // m = theta(3); //RF CG change for testing against TINSS - estimate M directly
	  log_avgrec = theta(4);
	  log_recinit = theta(5);
	  rho = theta(6);
    cout<<"rho="<<rho<<endl;
	  varphi = theta(7);

	  if(verbose)
      cout<<"**** Ok after initParameters ****"<<endl;
  }

FUNCTION dvar_vector cubic_spline(const dvar_vector& spline_coffs)
  {
	RETURN_ARRAYS_INCREMENT();
	int nodes=size_count(spline_coffs);
	dvector ia(1,nodes);
	dvector fa(sage,nage);
	ia.fill_seqadd(0,1./(nodes-1));
	fa.fill_seqadd(0,1./(nage-sage));
	vcubic_spline_function ffa(ia,spline_coffs);
	RETURN_ARRAYS_DECREMENT();

	//some testing here
	/*dvar_vector spline_nodes(1,nodes);
		spline_nodes.fill_seqadd(-0.5,1./(nodes-1));
		cout<<spline_nodes<<endl;
		vcubic_spline_function test_ffa(ia,spline_nodes);
		cout<<test_ffa(fa)<<endl;
		exit(1);
  */
	  return(ffa(fa));
  }

FUNCTION dvar_matrix cubic_spline_matrix(const dvar_matrix& spline_coffs)
  {
	  RETURN_ARRAYS_INCREMENT();
	  int nodes= spline_coffs.colmax()-spline_coffs.colmin()+1;
	  int rmin = spline_coffs.rowmin();
	  int rmax = spline_coffs.rowmax();

	  dvector ia(1,nodes);
	  dvector fa(sage,nage);
	  ia.fill_seqadd(0,1./(nodes-1));
	  //fa.fill_seqadd(sage,1);
	  fa.fill_seqadd(0,1./(nage-sage));
	  vcubic_spline_function_array fna(rmin,rmax,ia,spline_coffs);
	  RETURN_ARRAYS_DECREMENT();
	  return(fna(fa));
  }

FUNCTION calcSelectivities
    {
  	/*
  		This function loops over each ngear and calculates the corresponding
  		selectivity for that gear type. It first uses a switch statement 
  		to calculate selectivity curves based on isel_type where:
  		1) logistic selectivity with 2 parameters
  		2) age-specific selectivity coefficients with (nage-sage) parameters
  		   and the last two age-classes are assumed to have the same selectivity.
  		3) a reduced age-specific parameter set based on a bicubic spline.
  		4) Time varying cubic spline.
  		5) Time varying bicubic spline (2d version)
  		6) Fixed logistic
  		7) logistic selectivity based on relative changes in mean weight at age
  		
  		Following the initializatoin of the selectivity curves, time-varying 
  		considerations are implemented.
  		
  		TODO: Add penality (10.*square(avg_log_sel)) to objective function 
  		in cases where estimating sel_coffs to mimic an init_bounded_dev vector.
  		
  		CHANGED: Problem with case 7: turns out to be a random walk, so there
  		is changes in selectivity even with constant growth.  Remove the
  		random walk outside the switch statement.
  	
  	*/
  	int i,j,k;
  	double tiny=1.e-10;
  	dvariable p1,p2,p3;
  	dvar_vector age_dev=age;
  	dvar_matrix t1;
  	dvar_matrix tmp(syr,nyr-1,sage,nage);
  	dvar_matrix tmp2(syr,nyr,sage,nage);
  	dvar_matrix ttmp2(sage,nage,syr,nyr);
  	//jlog_sel.initialize();
  	log_sel.initialize();
  	avg_log_sel.initialize();
  	
  	for(j=1;j<=ngear;j++)
  	{
  		tmp.initialize(); tmp2.initialize();
  		dvector iy(1,yr_nodes(j));
  		dvector ia(1,age_nodes(j));
  		
  		switch(isel_type(j))
  		{
  			case 1:
  				// logistic selectivity for case 1 or 6
  				p1 = mfexp(sel_par(j,1,1));
  				p2 = mfexp(sel_par(j,1,2));
  				for(i=syr; i<=nyr; i++)
  				{
  					log_sel(j)(i) = log( plogis(age,p1,p2)+tiny );
  				}
  				break;
  			
  			case 6:
  				// logistic selectivity for case 1 or 6
  				p1 = mfexp(sel_par(j,1,1));
  				p2 = mfexp(sel_par(j,1,2));
  				for(i=syr; i<=nyr; i++)
  				{
  					log_sel(j)(i) = log( plogis(age,p1,p2) );
  				}
  				break;
  				
  			case 2:		
  				// age-specific selectivity coefficients
  				if(verbose) cout<<"age-specific sel"<<endl;
  				for(i=syr; i<=nyr; i++)
  				{
  					for(k=sage;k<=nage-1;k++)
  					log_sel(j)(i)(k) = sel_par(j)(1)(k-sage+1);
  					log_sel(j)(i,nage) = log_sel(j)(i,nage-1);
  				}
  				break;
  				
  			case 3:		
  				// cubic spline
  				log_sel(j)(syr)=cubic_spline( sel_par(j)(1) );
  				for(i=syr; i<nyr; i++)
  				{
  					log_sel(j)(i+1) = log_sel(j)(i);
  				}
  				break;
  				
  			case 4:		
  				// time-varying cubic spline every year
  				for(i=syr; i<=nyr; i++)
  				{
  					log_sel(j)(i) = cubic_spline(sel_par(j)(i-syr+1));
  				}
  				
  				//jlog_sel(j) = cubic_spline(sel_par(j)(1));
  				//t1 = cubic_spline_matrix(sel_par(j).sub(2,jsel_npar(j)));
  				//for(i = t1.indexmin(); i <= t1.indexmax(); i++)
  				//{
  				//	tmp( syr+(i-t1.indexmin()) ) = t1(i);
  				//}
  				break;
  				
  			case 5:		
  				// time-varying bicubic spline
  				ia.fill_seqadd( 0,1./(age_nodes(j)-1) );
  				iy.fill_seqadd( 0,1./(yr_nodes(j)-1) );
  				// bicubic_spline function is located in stats.cxx library
  				bicubic_spline( iy,ia,sel_par(j),tmp2 );
  				log_sel(j) = tmp2; 
  				break;
  				
  			case 7:
  				// time-varying selectivity based on deviations in weight-at-age
  				// CHANGED This is not working and should not be used. (May 5, 2011)
  				// SJDM:  I was not able to get this to run very well.
  				// AUG 5, CHANGED so it no longer has the random walk component.
  				p1 = mfexp(sel_par(j,1,1));
  				p2 = mfexp(sel_par(j,1,2));
  				
  				for(i = syr; i<=nyr; i++)
  				{
  					dvar_vector tmpwt=log(wt_obs(i)*1000)/mean(log(wt_obs*1000.));
  					log_sel(j)(i) = log( plogis(tmpwt,p1,p2)+tiny );
  				}	 
  				break;
  				
  			case 8:
  				//Alternative time-varying selectivity based on weight deviations (wt_dev)
  				//wt_dev is a matrix(syr,nyr+1,sage,nage)
  				//p3 is the coefficient that describes variation in log_sel.
  				p1 = mfexp(sel_par(j,1,1));
  				p2 = mfexp(sel_par(j,1,2));
  				p3 = sel_par(j,1,3);
  				
  				for(i=syr; i<=nyr; i++)
  				{
  					tmp2(i) = p3*wt_dev(i);
  					log_sel(j)(i) = log( plogis(age,p1,p2)+ tiny ) + tmp2(i);
  				}
  				break;
  				
  			default:
  				log_sel(j)=0;
  				break;
  				
		} 
		//log_sel(j)(syr) = jlog_sel(j)+tmp2(syr);
		//for(i=syr;i<nyr;i++)
		//{
		//	log_sel(j)(i+1)(sage,nage) = log_sel(j)(i)(sage,nage)+tmp(i)+tmp2(i+1);
		//}
		//subtract mean to ensure mean(exp(log_sel))==1
		//substract max to ensure exp(log_sel) ranges from 0-1
		
		//Correction not applied for cases 1,6,7 - was returning values greater than one when exponentiated
		if(isel_type(j)==2 || isel_type(j)==3|| isel_type(j)==4|| isel_type(j)==5) {
		//subtract mean to ensure sum(log_sel)==0
			for(int i=syr;i<=nyr;i++)
				log_sel(j)(i) -= log( mean(mfexp(log_sel(j)(i)))+tiny );
				//log_sel(j)(i) -= log(mean(mfexp(log_sel(j)(i))));
				//log_sel(j)(i) -= log(max(mfexp(log_sel(j)(i))));
		}
		//cout<<"log_sel"<<log_sel<<endl;
		//cout<<"log_sel \t"<<j<<"\n"<<log_sel(j)<<"\n \n"<<endl;
		//testing bicubic spline  (SM Checked OCT 25,2010.  Works on the example below.)
		/*ia.fill_seqadd(0,1./(age_nodes(j)-1));
				iy.fill_seqadd(0,1./(yr_nodes(j)-1));
				dvar_matrix tn(1,age_nodes(j),1,yr_nodes(j));
				tn.colfill_seqadd(1,-.5,0.1);
				tn.colfill_seqadd(2,-.4,0.1);
				bicubic_spline(iy,ia,tn,tmp2);
				cout<<ia<<endl;
				cout<<iy<<endl;
				cout<<tn<<endl;
				cout<<tmp2<<endl;
				exit(1);*/
	}
	
	if(verbose)cout<<"**** Ok after calcSelectivities ****"<<endl;

  }

FUNCTION calcTotalMortality
  {
	/*
	This routine calculates fishing mortality, total mortality
	and annaul survival rates (exp(-Z)) for each age in each
	year.
	There is a complication in that if there is a survey gear
	then the fishing mortality rate is set to an extremely small
	value: exp(-70.)~3.975e-31, which is effectively 0.
	The above issue is no longer an issue b/c now estimating Fs
	for years where there is non-zero catch.
	SJDM.  Dec 24, 2010.  Adding time-varying natural mortality
	CHANGED May 20, 2011  Add cubic spline to the time-varying natural mortality
	*/
	int i,k,ki;
	dvariable ftmp;
	F.initialize();
	ft.initialize();
	log_ft.initialize();

	//Fishing mortality
	ki=1;
	for(k=1;k<=ngear;k++)
	{
		for(i=syr;i<=nyr;i++)
		//for(i=f_syr(k);i<=f_nyr(k);i++)
		{
			/*if(i==0) break;
						if(active(log_avg_f(k)))
							log_ft(k,i)=log_avg_f(k)+log_ft_devs(k,i);
						else log_ft(k,i)=-70.;
      */
			ftmp=0;
			if(obs_ct(k,i)>0)
				ftmp = mfexp(log_ft_pars(ki++));
			ft(k,i)=ftmp;
			//F(i)+=mfexp(log_ft(k,i)+log_sel(k)(i));
			F(i)+=ftmp*mfexp(log_sel(k)(i));
		}
	}
	//Natural mortality (year and age specific)
	//M_tot(syr,nyr,sage,nage);
	M_tot = m;
	// Cubic spline to interpolate log_m_devs (log_m_nodes)
	log_m_devs = 0.;
	if(active(log_m_nodes))
	{
		int nodes = size_count(log_m_nodes);
		dvector im(1,nodes);
		dvector fm(syr+1,nyr);
		im.fill_seqadd(0,1./(nodes-1));
		fm.fill_seqadd(0,1./(nyr-syr));
		vcubic_spline_function m_spline(im,log_m_nodes);
		//m_spline(fm);
		log_m_devs = m_spline(fm);
	}
	// Random walk in natural mortality.
	for(i=syr;i<=nyr;i++)
	{
		// if(active(log_m_devs)&&i>syr)
		if(active(log_m_nodes)&&i>syr)
		{
			M_tot(i)=M_tot(i-1)*exp(log_m_devs(i));
		}
	}
	m_bar = mean(M_tot);
	Z=M_tot+F;
	S=mfexp(-Z);
	if(verbose) cout<<"**** OK after calcTotalMortality ****"<<endl;
  }

FUNCTION calcNumbersAtAge
  {
	/*
		TODO Need to check the difference between the initialization 
		of the numbers at age here at the margins in comparison to the
		simulation model.
	*/
	int i,j;
	N.initialize();
	/*
	//RF BACK THIS OUT!
	//RF temp change for testing tinss
	for(i=sage+1; i<=nage;i++){
		 init_log_rec_devs(i)=log_rec_devs(syr-nage+i);
	}
	//end temp change
	*/
	if(cntrl(5)){	//If initializing in at unfished conditions
		log_rt(syr) = log(ro);
		for(j=sage;j<=nage;j++)
		{
			N(syr,j)=ro*exp(-m_bar*(j-1.));
		}
	}
	else{			//If starting at unfished conditions
		log_rt(syr) = log_avgrec+log_rec_devs(syr);
		N(syr,sage)=mfexp(log_rt(syr));
		for(j=sage+1;j<=nage;j++)
		{
			N(syr,j)=mfexp(log_recinit+init_log_rec_devs(j))*exp(-m_bar*(j-sage));
		}
	}
	N(syr,nage)/=(1.-exp(-m_bar));
	//initial number of sage recruits from year syr+1, nyr;
	for(i=syr+1;i<=nyr;i++){
		log_rt(i)=log_avgrec+log_rec_devs(i);
		N(i,sage)=mfexp(log_rt(i));
	}
	N(nyr+1,sage)=mfexp(log_avgrec);

	/*
	for(j=sage;j<=nage;j++)
	{
		if(cntrl(5))  //if starting at unfished state
		{
			N(syr,j)=ro*exp(-m_bar*(j-1));
		}
		else{
			log_rt(syr-j+sage)=log_avgrec+log_rec_devs(syr-j+sage);
			N(syr,j)=mfexp(log_rt(syr-j+sage))*exp(-m_bar*(j-sage));
		}
	}*/
	for(i=syr;i<=nyr;i++)
	{
		N(i+1)(sage+1,nage)=++elem_prod(N(i)(sage,nage-1),S(i)(sage,nage-1))+1.e-10;
		N(i+1,nage)+=N(i,nage)*S(i,nage);
	}
	if(verbose)cout<<"**** Ok after calcNumbersAtAge ****"<<endl;
  }

FUNCTION calcAgeProportions
  {
	/*This function loops over each gear and year
	and calculates the predicted proportions at age
	sampled based on the selectivity of that gear and
	the numbers-at-age in the population.*/
	int i,k,iyr,ig;
	for(k=1;k<=na_gears;k++)
	{
		for(i=1;i<=na_nobs(k);i++)
		{
			iyr=A(k,i,a_sage(k)-2);	//index for year
			ig=A(k,i,a_sage(k)-1);	//index for gear
			if(iyr>nyr)break;		//trap for retrospective analysis
			A_nu(k,i,a_sage(k)-2)=iyr;
			A_nu(k,i,a_sage(k)-1)=ig;
			Ahat(k,i,a_sage(k)-2)=iyr;
			Ahat(k,i,a_sage(k)-1)=ig;
			Ahat(k)(i)(a_sage(k),a_nage(k))=Chat(k)(iyr)(a_sage(k),a_nage(k))
										/sum(Chat(k)(iyr)(a_sage(k),a_nage(k)));
		}
	}
	if(verbose)cout<<"**** Ok after calcAgeProportions ****"<<endl;

  }

FUNCTION calcFisheryObservations
  {
	/*
	Dec 6, 2010.  Modified ct calculations to include
				  empirical weight at age data (wt_obs);
	Jan 16, 2011. 	modified this code to get age-comps from surveys, rather than 
					computing the age-comps in calc_fisheries_observations
	*/
	/*
		FIXED Reconcile the difference between the predicted catch here and in the simulation model.
	*/
	int i,k;
	ct.initialize();
	for(i=syr;i<=nyr;i++)
	{
		for(k=1;k<=ngear;k++)
		{
			dvar_vector log_va=log_sel(k)(i);// - log(mean(mfexp(log_sel(k)(i))));
			//dvar_vector fa=mfexp(log_ft(k,i)+log_va);
			//SJDM Jan 16, 2011 Modification as noted above.
			dvar_vector fa=ft(k,i)*mfexp(log_va);
			if(obs_ct(k,i)>0)
			{/*If there is a commercial fishery, then calculate the
			   catch-at-age (in numbers) and total catch (in weight)*/
				Chat(k,i)=elem_prod(elem_prod(elem_div(fa,Z(i)),1.-S(i)),N(i));//+1.e-10;
				ct(k,i) = Chat(k,i)*wt_obs(i);
			}
			else
			{/*If there is no commercial fishery the set Chat equal to
			   the expected proportions at age.*/
				fa = mfexp(log_va);
				Chat(k,i)=elem_prod(elem_prod(elem_div(fa,Z(i)),1.-S(i)),N(i));
			}
			/*
			Changed this on Jan 16, 2011 as it was preventing
			convergence for simuation with zero error due to the tiny
			constant added to F.
			need to add a tiny constant to deal with catch-age data
			for non-extractive survey age comps
			dvar_vector fa=ft(k,i)*mfexp(log_va)+1.e-30;
			//Catch-at-age by commercial gear
			if(fsh_flag(k))
				Chat(k,i)=elem_prod(elem_prod(elem_div(fa,Z(i)),1.-S(i)),N(i));
			//Catch weight by gear
			//ct(k,i)=Chat(k,i)*wa; 
			ct(k,i)=Chat(k,i)*wt_obs(i);  //SM Dec 6, 2010*/
		}
	}
	if(verbose)cout<<"**** Ok after calcFisheryObservations ****"<<endl;

  }

FUNCTION calc_survey_observations
  {
	/*This code needs to be modified to accomodate
	multiple surveys or block changes in survey q.
	Oct 31, 2010, added retrospective counter.
	Nov 22, 2010, adding multiple surveys. Still need to check with retrospective option
	Nov 30, 2010, adjust the suvery biomass by the fraction of Z that has occurred 
	when the survey was conducted. For herring spawning biomass this would be after the 
	fishery has taken place.
	Dec 6, 2010, modified predicted survey biomass to accomodate empirical weight-at-age 
	data (wt_obs).
	May 11, 2011.  CHANGED Vivian Haist pointed out an error in survey biomass comparison.
	The spawning biomass was not properly calculated in this routine. I.e. its different 
	than the spawning biomass in the stock-recruitment routine. (Based on fecundity which
	changes with time when given empirical weight-at-age data.)
	*/
	/*
		CHANGED add capability to accomodate priors for survey q's.
		DONE
	*/
	int i,j,ii,k;
	//survey abudance index residuals
	epsilon.initialize();
	pit.initialize();
	for(i=1;i<=nit;i++)
	{
		int nx=0;		//counter for retrospective analysis
		dvar_matrix V(1,nit_nobs(i),sage,nage);  //vulnerable units for survey comparison
		V.initialize();
		//dvar_matrix VB(1,nit_nobs(i),sage,nage); //vulnerable biomass
		for(j=1;j<=nit_nobs(i);j++)
		{
			ii=iyr(i,j);
			k=igr(i,j);
			if(ii>nyr) break;	//trap for retrospective analysis.
			dvar_vector log_va=log_sel(k)(ii);
			//Adjust survey biomass by the fraction of the mortality that
			//occurred during the time of the survey.
			dvar_vector Np = elem_prod(N(ii),exp( -Z(ii)*it_timing(i,j) ));
			//V(j)=elem_prod(N(ii),mfexp(log_va));
			switch(survey_type(i))
			{
				case 1:
					V(j)=elem_prod(Np,mfexp(log_va));
				break;
				case 2:
					V(j)=elem_prod(elem_prod(Np,mfexp(log_va)),wt_obs(ii));
				break;
				case 3:
					V(j)=elem_prod(Np,fec(ii));
				break;
			}
			//VB(j)=elem_prod(V(j),wt_obs(ii));		//SM Dec 6, 2010
			//If the survey is a spawn index, then need to account
			//for changes in fecundity.
			if(iyr(i,j)<=nyr) nx++;
		}
		dvar_vector t1 = rowsum(V);//V*wa;
		//See Ludwig & Walters 1994
		//Note this is incorrect if each survey has different weights.
		dvar_vector zt=log(it(i).sub(1,nx))-log(t1(1,nx));
		//cout<<"zt\n"<<t1(1,nx)<<endl;
		epsilon(i).sub(1,nx) = zt-mean(zt);
		q(i) = exp(mean(zt));
		pit(i).sub(1,nx)=t1(1,nx)*q(i);	//predicted index
		//TODO, this might not be working correctly, simulation test it.
		if(q_prior(i)==2)
		{
			//random walk in q
			epsilon(i)=0;
			dvar_vector fd_zt=first_difference(zt);
			epsilon(i).sub(1,nx-1) = fd_zt-mean(fd_zt);
			//dvar_vector qt(1,nx);
			qt(i,1) = exp(zt(1));
			for(j=2;j<=nx;j++)
				qt(i,j) = qt(i,j-1)*exp(fd_zt(j-1));
			pit(i).sub(1,nx)=elem_prod(t1(1,nx),qt(i)(1,nx));
		}
	}
	if(verbose)cout<<"**** Ok after calc_survey_observations ****"<<endl;
  }

FUNCTION calc_stock_recruitment
  {
	/*
	The following code is used to derive unfished
	spawning stock biomass bo and the stock-
	recruitment parameters for the:
	Beverton-Holt Model
		-Rt=k*Ro*St/(Bo+(k-1)*St)*exp(delta-0.5*tau*tau)
	Ricker Model
		-Rt=so*St*exp(-beta*St)*exp(delta-0.5*tau*tau)
	Dec 6, 2010.  Modified code to allow empirical weight-at-age data
				This required a fecundity-at-age matrix.  Need to
				project spawning biomass into the future.
	CHANGED bo should be a function of the average natural mortality
	TODO update phib calculation if age-specific M's are used.
	May 6, 2010.  Changed phib calculation based on average M
	in cases where M varies over time. Otherwise bo is biased based
	on the initial value of M.
	CHANGED Need to adjust spawning biomass to post fishery numbers.
	CHANGED Need to adjust spawners per recruit (phib) to average fecundity.
	*/
	int i;
	dvariable tau = (1.-rho)/varphi;
	dvar_vector tmp_rt(syr+sage,nyr);
	dvar_vector lx(sage,nage); lx=1;
	for(i=sage+1;i<=nage;i++) lx(i)=lx(i-1)*exp(-m_bar);
	lx(nage)/=(1.-exp(-m_bar));
	dvariable phib = (lx*exp(-m_bar*cntrl(13))) * avg_fec;	//SM Dec 6, 2010
	// No redeclaration of so and beta - they are declared in parameter section - makes them available to REPORT section
	so = kappa/phib;		//max recruits per spawner
	bo = ro*phib;  			//unfished spawning biomass
	if(verbose) cout<<"Init bo = "<<bo<<endl;
	//sbt=rowsum(elem_prod(N,fec));			//SM Dec 6, 2010
	//CHANGED adjusted spawning biomass downward by ctrl(12)
	for(i=syr;i<=nyr;i++)
	{
		sbt(i) = elem_prod(N(i),exp(-Z(i)*cntrl(13)))*fec(i);
	}
	sbt(nyr+1) = N(nyr+1)*fec(nyr+1);
	dvar_vector tmp_st=sbt(syr,nyr-sage).shift(syr+sage);
	switch(int(cntrl(2)))
	{
		case 1:
			//Beverton-Holt model
			beta = (kappa-1.)/bo;
			tmp_rt=elem_div(so*tmp_st,1.+beta*tmp_st);
			break;
		case 2:
			//Ricker model
			//dvariable so=kappa/phib;
			beta=log(kappa)/bo;
			tmp_rt=elem_prod(so*tmp_st,exp(-beta*tmp_st));
		break;
	}
	//residuals in stock-recruitment curve
	rt=exp(log_rt(syr+sage,nyr));//trans(N)(1)(syr+1,nyr);
	delta = log(rt)-log(tmp_rt)+0.5*tau*tau;
	if(verbose)cout<<"**** Ok after calc_stock_recruitment ****"<<endl;
  }

FUNCTION calc_objective_function
  {
	//Dec 20, 2010.  SJDM added prior to survey qs.
	/*q_prior is an ivector with current options of 0 & 1.
	0 is a uniform density (ignored) and 1 is a normal
	prior density applied to log(q).*/

	/*
	There are several components to the objective function
	Likelihoods:
		-1) likelihood of the catch data
		-2) likelihood of the survey abundance index
		-3) likelihood of the survey age comps
		-4) likelihood of the fishery age comps
		-5) likelihood for stock-recruitment relationship
		-6) likelihood for fishery selectivities
	*/
	int i,j,k;
	double o=1.e-10;

	dvar_vector lvec(1,6); lvec.initialize();
	nlvec.initialize();
	//1) likelihood of the catch data (retro)
	double sig_c =cntrl(3);
	if(last_phase())sig_c=cntrl(4);
	if(active(log_ft_pars))
	for(k=1;k<=ngear;k++){
		for(i=syr;i<=nyr;i++)
		{
			if(obs_ct(k,i)!=0)
				nlvec(1,k)+=dnorm(log(ct(k,i)),log(obs_ct(k,i)),sig_c);
		}
		//if(active(log_ft_pars))
		//	nlvec(1,k)=dnorm(log(obs_ct(k).sub(syr,nyr)+o)-log(ct(k).sub(syr,nyr)+o),sig_c);
	}
	//2) likelihood of the survey abundance index (retro)
	for(k=1;k<=nit;k++)
	{
		dvar_vector sig = (rho/varphi)/it_wt(k);
		nlvec(2,k)=dnorm(epsilon(k),sig);
	}
	//3) likelihood for age-composition data
	for(k=1;k<=na_gears;k++)
	{
		if(na_nobs(k)>0){
			int naa=0;
			int iyr;
			//retrospective counter
			for(i=1;i<=na_nobs(k);i++)
			{
				iyr=A(k,i,a_sage(k)-2);	//index for year
				if(iyr<=nyr)naa++;
			}
			dmatrix O=trans(trans(A(k)).sub(a_sage(k),a_nage(k))).sub(1,naa);
			dvar_matrix P=trans(trans(Ahat(k)).sub(a_sage(k),a_nage(k))).sub(1,naa);
			//dvar_matrix nu=trans(trans(Ahat(k)).sub(a_sage(k),a_nage(k))).sub(1,naa); //residuals
			dvar_matrix nu(O.rowmin(),O.rowmax(),O.colmin(),O.colmax()); //residuals
			nu.initialize();
			//CHANGED add a switch statement here to choose form of the likelihood
			//RF Added 3rd switch: 3=no likelihood (for "ageless" model) - currently still need a dummy matrix to avoid read data errors
			switch(int(cntrl(14)))
			{
				case 1:
					nlvec(3,k) = dmvlogistic(O,P,nu,age_tau2(k),cntrl(6));
				break;
				case 2:
					nlvec(3,k) = dmultinom(O,P,nu,age_tau2(k),cntrl(6));
				break;
				case 3:
					nlvec(3,k) = 0.; //rf added this July 23 2012 for 'ageless' model
				break;
			}
			//if(last_phase) cout<<nlvec(3,1)<<endl;
			for(i=1;i<=naa/*na_nobs(k)*/;i++)
			{
				iyr=A(k,i,a_sage(k)-2);	//index for year
				A_nu(k)(i)(a_sage(k),a_nage(k))=nu(i);
			}
		}
	}
	//4) likelihood for stock-recruitment relationship
	dvariable tau = (1.-rho)/varphi;
	if(active(theta(1)))
		nlvec(4,1)=dnorm(delta,tau);
	//5-6) likelihood for selectivity paramters
	for(k=1;k<=ngear;k++)
	{
		if(active(sel_par(k))){
			//if not using logistic selectivity then
			//CHANGED from || to &&  May 18, 2011 Vivian
			if( isel_type(k)!=1 && isel_type(k)!=7 && isel_type(k)!=8 )
			{
				for(i=syr;i<=nyr;i++)
				{
					//curvature in selectivity parameters
					dvar_vector df2=first_difference(first_difference(log_sel(k)(i)));
					nlvec(5,k)+=sel_2nd_diff_wt(k)/(nage-sage+1)*norm2(df2);
					//penalty for dome-shapeness
					for(j=sage;j<=nage-1;j++)
						if(log_sel(k,i,j)>log_sel(k,i,j+1))
							nlvec(6,k)+=sel_dome_wt(k)*square(log_sel(k,i,j)-log_sel(k,i,j+1));
				}
			}
		}
	}
	// CONSTRAINT FOR SELECTIVITY DEV VECTORS
	// Ensure vector of sel_par sums to 0. (i.e., a dev_vector)
	// TODO for isel_type==2 ensure mean 0 as well (ie. a dev_vector)
	for(k=1;k<=ngear;k++)
	{
		if( active(sel_par(k)) && isel_type(k)!=1 && isel_type(k)!=7 && isel_type(k)!=8 )
		{
			dvariable s=0;
			if(isel_type(k)==5)  //bicubic spline version ensure column mean = 0
			{
				dvar_matrix tmp = trans(sel_par(k));
				for(j=1;j<=tmp.rowmax();j++)
				{
					s=mean(tmp(j));
					lvec(1)+=1000.*s*s;
				}
			}
			if(isel_type(k)==4 || isel_type(k)==3)
			{
				dvar_matrix tmp = sel_par(k);
				for(j=1;j<=tmp.rowmax();j++)
				{
					s=mean(tmp(j));
					lvec(1)+=1000.*s*s;
				}
			}
		}
	}
	/*if(active(spline_coffs)||active(spline2_coffs))
		{
			//lvec(6)=1.*norm2(spline2_coffs);
			for(i=syr;i<=nyr;i++)
			{
				//curvature in selectivity parameters
				dvar_vector df2=first_difference(first_difference(log_sel(i)));
				lvec(6)+=50.0/(nyr-syr+1)*norm2(df2);
				//penalty for dome-shapeness
				for(j=sage;j<=nage-2;j++)
					if(log_sel(i,j)>log_sel(i,j+1))
						lvec(6)+=3.125*square(log_sel(i,j)-log_sel(i,j+1));
			}
			//penalty for random walk in age-changes
			for(j=sage;j<=nage;j++){
				lvec(6)+=0.*norm2(first_difference(trans(log_sel)(j)));
			}
		}
		else
		{
			dvar_vector df2=first_difference(first_difference(log_sel(syr)));
			lvec(6)=50./(nyr-syr+1)*norm2(df2);
		}*/
	/*
	@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
	~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	 PRIORS for estimated model parameters from the control file
	~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
	*/
	//pvec(1)=dlnorm(msy,log(msyMean),msySd);
	dvariable ptmp; priors.initialize();
	for(i=1;i<=npar;i++)
	{
		if(active(theta(i)))
		{
			switch(theta_prior(i))
			{
			case 1:		//normal
				ptmp=dnorm(theta(i),theta_control(i,6),theta_control(i,7));
				break;
			case 2:		//lognormal CHANGED RF found an error in dlnorm prior. rev 116
				ptmp=dlnorm(theta(i),theta_control(i,6),theta_control(i,7));
				break;
			case 3:		//beta distribution (0-1 scale)
				double lb,ub;
				lb=theta_lb(i);
				ub=theta_ub(i);
				ptmp=dbeta((theta(i)-lb)/(ub-lb),theta_control(i,6),theta_control(i,7));
				break;
			case 4:		//gamma distribution
				ptmp=dgamma(theta(i),theta_control(i,6),theta_control(i,7));
				break;
			default:	//uniform density
				ptmp=1./(theta_control(i,3)-theta_control(i,2));
				break;
			}
			priors(i)=ptmp;
		}
	}
	// added prior for ghat (same prior for all gears)
	for(k=1;k<=ngear;k++)
	{
		if( isel_type(k)==1 || isel_type(k)==6 )
		{
			if(ghat_pswitch(k)) priors(npar+1)=dgamma(mfexp(sel_par(k,1,2)),ghat_p1,ghat_p2);
		}
	}
	
	//Priors for suvey q based on control file.
	qvec.initialize();
	for(i=1;i<=nits;i++)
	{
		if(q_prior(i)==1)
		{
			qvec(i)=dnorm(log(q(i)),q_mu(i),q_sd(i));
		}
	}
	/*
	The following are penalties that are invoked in early
	phases to ensure reasonable solutions are obtained,
	and to reduce the sensitivity of initial starting
	conditions.  Penalties include:
		-1) keep average fishing mortality rate near
			0.2 and in the last phase relax this constraint.
		-3) normal prior for log rec devs with std=50.
	*/
	pvec.initialize();
	//Penalties to regularize the solution for fishing mortality rates
	dvariable log_fbar = mean(log_ft_pars);
	if(last_phase())
	{
		pvec(1) = dnorm(log_fbar,log(cntrl(7)),cntrl(9));
		//Penalty for log_rec_devs (large variance here)
		pvec(4) = dnorm(log_rec_devs,2.0); //2.0
		pvec(5) = dnorm(init_log_rec_devs,2.0); //2.0
		// set same as TINSS for testing
		//pvec(4) = 0.01*norm2(log_rec_devs);
		//pvec(5) = 0.01*norm2(init_log_rec_devs);
	}
	else
	{
		pvec(1) = dnorm(log_fbar,log(cntrl(7)),cntrl(8));
		//Penalty for log_rec_devs (CV ~ 0.0707) in early phases
		pvec(4)=100.*norm2(log_rec_devs);//
		pvec(5)=100.*norm2(init_log_rec_devs); //
	}
	//Priors for deviations in natural mortality rates
	//if(active(log_m_devs))
	if(active(log_m_nodes))
	{
		double std_mdev = cntrl(11);
		dvar_vector fd_mdevs=first_difference(log_m_devs);
		pvec(2) = dnorm(fd_mdevs,std_mdev);
		pvec(2) += 0.5*norm2(log_m_nodes);
	}
	if(verbose)
	{
		cout<<"nlvec\t"<<nlvec<<endl;
		cout<<"lvec\t"<<lvec<<endl;
		cout<<"priors\t"<<priors<<endl;
		cout<<"penalties\t"<<pvec<<endl;
	}
	// add pf for priors
	Pf=sum(priors)+sum(qvec);
	Lf=sum(nlvec)+sum(lvec)+sum(pvec);
	//f=sum(nlvec)+sum(lvec)+sum(priors)+sum(pvec)+sum(qvec);
	f=Lf+Pf;
	nf++;
	if(verbose)cout<<"**** Ok after calc_objective_function ****"<<endl;
  }

FUNCTION void equilibrium(const double& fe,const double& ro, const double& kap, const double& m, const dvector& age, const dvector& wa, const dvector& fa, const dvector& va,double& re,double& ye,double& be,double& phiq,double& dphiq_df, double& dre_df, double& dspr_df,double& dsb_df, double& spre, double& ue)
  {
	/*
	This is the equilibrium age-structured model that is
	conditioned on fe (the steady state fishing mortality rate).
	In the case of multiple fisheries, fe is to be considered as the
	total fishing mortality rate and each fleet is given a specified
	allocation based on its selectivity curve.  The allocation to
	each fleet must be specified a priori.
	args:
	fe	-steady state fishing mortality
	ro	-unfished sage recruits
	kap	-recruitment compensation ration
	m	-instantaneous natural mortality rate
	age	-vector of ages
	wa	-mean weight at age
	fa	-mean fecundity at age
	va	-mean vulnerablity at age for fe gear.

	Modified args:
	re	-steady state recruitment
	ye	-steady state yield
	be	-steady state spawning biomass
	phiq		-per recruit yield
	dre_df		-partial of recruitment wrt fe
	dphiq_df	-partial of per recruit yield wrt fe

	FIXME add Ricker model to reference points calculations.
	FIXME partial derivatives for dphif_df need to be fixed when cntrl(12)>0.
	*/
	int i;
	int nage=max(age);
	int sage=min(age);
	dvector lx=pow(exp(-m),age-double(sage));
	lx(nage)/=(1.-exp(-m));
	dvector lz=lx;
	dvector za=m+fe*va;
	dvector sa=1.-exp(-za);
	dvector qa=elem_prod(elem_div(va,za),sa);
	double phie = lx*fa;		//eggs per recruit
	double so = kap/phie;
	double beta = (kap-1.)/(ro*phie);
	double dlz_df = 0, dphif_df = 0;
	dphiq_df=0; dre_df=0;
	for(i=sage; i<=nage; i++)
	{
		if(i>sage) lz[i]=lz[i-1]*exp(-za[i-1]);
		if(i>sage) dlz_df=dlz_df*exp(-za[i-1]) - lz[i-1]*va[i-1]*exp(-za[i-1]);
		if(i==nage){ //6/11/2007 added plus group.
					lz[i]/=(1.-mfexp(-za[i]));
					//dlz_df=dlz_df*mfexp(-za[i-1]) - lz[i-1]*va[i-1]*mfexp(-za[i-1])/(1.-mfexp(-za[i]))
					dlz_df=dlz_df/(1.-mfexp(-za[i]))
							-lz[i-1]*mfexp(-za[i-1])*va[i]*mfexp(-za[i])
					/((1.-mfexp(-za[i]))*(1.-mfexp(-za[i])));
				}
		dphif_df=dphif_df+fa(i)*dlz_df;
		dphiq_df=dphiq_df+wa(i)*qa(i)*dlz_df+(lz(i)*wa(i)*va(i)*va(i))/za(i)*(exp(-za[i])-sa(i)/za(i));
	}
	//CHANGED need to account for fraction of mortality that occurs
	//before the spawning season in the recruitment calculation.
	//cout<<"lz\t"<<elem_prod(lz,exp(-za*cntrl(13)))<<endl;
	//exit(1);
	//double phif = lz*fa;

	double phif = elem_prod(lz,exp(-za*cntrl(13)))*fa;
	phiq=sum(elem_prod(elem_prod(lz,wa),qa));
	re=ro*(kap-phie/phif)/(kap-1.);
	//cout<<fe<<" spr ="<<phif/phie<<endl;
	//spr_test=phif/phie;
	if(re<=0) re=0.;
	dre_df=(ro/(kap-1.))*phie/square(phif)*dphif_df;
	ye=fe*re*phiq;
	be=re*phif;	//spawning biomass

	double be3=0.;
	for(i=3;i<=nage;i++) be3 += re*lz(i)*wa(i);
	ue =ye/be3;
	spre = value(-ro/((kappa-1)*re-ro*kappa));
	dspr_df=dphif_df/phie;
	dsb_df=re*dphif_df+phif*dre_df;
  }

FUNCTION void calc_reference_points()
  {
	/**
	\file ccam.tpl
	\author Steven Martell
	\\Added to by Robyn Forrest
	Uses Newton_Raphson method to determine Fmsy and MSY
	based reference points.

	Code check: appears to find the correct value of MSY
	in terms of maximizing ye.  Check to ensure rec-devs
	need a bias correction term to get this right.

	Code check: In management mode, with verbose =1, the equilibrium fe and ye values
	tested to return the same values as estimated fmsy and msy - check of back calculation

	Added code to calculate spr based reference points for Pacific hake assessment 2012

	Modification for multiple fleets:
		Need to pass a weighted average vector of selectivities
		to the equilibrium routine, where the weights for each
		selectivity is based on the allocation to each fleet.

		Perhaps as a default, assign an equal allocation to each
		fleet.  Eventually,user must specify allocation in
		control file.
		Use selectivity in the terminal year to calculate reference
		points.
	*/
	int i,j;
	double re,ye,be,de,phiq,dphiq_df,dre_df,fe; //RF added de (depletion)
	double dye_df,ddye_df,spr;
	fe = value(m_bar); //initial
	//cout<<"Initial fe = "<<fe<<endl<<endl; 

	// add this to arguments of equilibrium function - for calculating F40
	double dspr_df;
	double dsb_df;
	double spr4; // another name for the spawning potential ratio - keep separate from other calcs tmp for getting f40

	/*Calculate average vulnerability (among gears) in final year*/
	dvector va_bar(sage,nage);
	va_bar.initialize();
	/*CHANGED user now specifies allocation for MSY based reference points 
	in the data file.  Used to be fsh_flag, but now is an allocation for gear k*/
	//dvector allocation(1,ngear);
	//allocation = dvector(fsh_flag/sum(fsh_flag));

	/*CHANGED Allow for user to specify allocation among gear types.*/
	for(j=1;j<=ngear;j++)
	{
		va_bar+=allocation(j)*value(exp(log_sel(j)(nyr)));
	}

	/*CHANGED Changed equilibrium calculations based on average m */
	// FIXME: change Bmsy calculations to be based on average fecundity & weight at age. 

	//1. Loop over test values of f to find f that maximises equilibrium yield (in management mode, this is just a test of the equilbrium code)
	//MSY reference points
		for(i=1;i<=25;i++)
		{
			// use average weight and fecundity to be consistet with calculation of bo as well as calculation of fmsy and msy in management oriented mode
			equilibrium(fe,value(ro),value(kappa),value(m_bar),age,avg_wt,avg_fec,va_bar,re,msy_y,msy_sb,phiq,dphiq_df,dre_df, dspr_df,dsb_df, msy_spr, msy_u);
			//equilibrium(fe,value(ro),value(kappa),value(m_bar),age,wt_obs(nyr), fec(nyr),va_bar,re,ye,be,phiq,dphiq_df,dre_df, dspr_df,dsb_df);
			dye_df = re*phiq+fe*phiq*dre_df+fe*re*dphiq_df;
			ddye_df = phiq*dre_df + re*dphiq_df;
			fe = fe - dye_df/ddye_df;
			if(verbose) cout<<"fe\t"<<fe<<"\t"<<dye_df<<"\t"<<ye<<endl;
			if(sfabs(dye_df)<1.e-7)break;
		}
		// TESTING MANAGEMENT-ORIENTED CODE
		if(useManagementParameters) {
			if(verbose){
				cout<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
				cout<<"Testing Management Code"<<endl;
				cout<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
				cout<<"Estimated management parameters"<<endl;
				cout<<"Fmsy = "<<fmsy<<"\t"<<"MSY = "<<msy<<endl;
				cout<<"Back-calculated management parameters"<<endl;
				cout<<"Fe = "<<fe<<"\t"<<"Ye = "<<ye<<endl<<endl;
				if(sfabs(fmsy-fe) > 1.e-4) cout<<"~~~~~~~~~~FMSY MISMATCH!!~~~~~~~~~~~"<<endl;
				if(sfabs(msy-ye) > 1.e-4) cout<<"~~~~~~~~~~MSY MISMATCH!!~~~~~~~~~~~"<<endl<<endl;
			}
		}

	if(!useManagementParameters) {
		fmsy=fe;
		msy=msy_y;
		
	}
	//cout<<"msy "<<msy<<endl; exit(1);
	//2. Run equilibrium again with fmsy. This finds bmsy and also spr at fmsy
	//value(fmsy) because fmsy has been changed to parameter. Also use average weight and fecundity to be consistent with calculation of bo as well as calculation of fmsy and msy in management oriented mode
	equilibrium(value(fmsy),value(ro),value(kappa),value(m_bar),age,avg_wt,avg_fec,va_bar,re,ye,be,phiq,dphiq_df,dre_df,dspr_df, dsb_df,msy_spr, msy_u);
	bmsy=be;
	//spr_msy = value(-ro/((kappa-1)*re-ro*kappa)); //this line is now redundant as it is passed by reference above (as msy_spr) - TESTED RF
	/*TODO print this to the REPORT file for plotting.*/
	/*SM Loop over discrete value of fe and ensure above code is
	finding the correct value of msy.*/

	//3. Run the equilibrium code over a range of f and keep values - tests that the value of fmsy found above really does maximise yield - keeps track of spr and eqm recruitment
	
	if(!mceval_phase()){
		ofstream report_file("ccam.eql");

		if(report_file.is_open())
		{
			report_file<<"index\t fe \t ye \t be \t de \t re \t spr\n"; //RF added de (depletion)

				fe = 0; i=0;
			double spre, ue;
			while(i < 1500)
			{
				// use average weight and fecundity to be consistent with calculation of bo as well as calculation of fmsy and msy in management oriented mode
				equilibrium(fe,value(ro),value(kappa),value(m_bar),age,avg_wt,avg_fec,va_bar,re,ye,be,phiq,dphiq_df,dre_df,dspr_df, dsb_df, spre, ue);
				//End change
				if(re<=0)break;
				if(verbose){
					if(i==0) {
						cout<<"Ref bo = "<<bo<<endl;
						cout<<"Ref be = "<<be<<endl<<endl;}}//End if end if
				spr = value(-ro/((kappa-1)*re-ro*kappa));
				de = be/value(bo);
				report_file<<i++<<"\t"<<fe<<"\t"<<ye<<"\t"<<be<<"\t"<<de<<"\t";
				report_file<<re<<"\t"<<spr<<endl;
				fe += 0.001;
			}
		}
	}
	//4. Run the equilibrium code again under the historical estimated ft to get spr for each year
	//For Pacific hake, ignore second row (survey)
	//For multiple fleets, need to be smarter than this bc need to consider f over whole fishery
	if(last_phase()){
		double Ft;
		double spre, ue;
		for(i=syr;i<=nyr;i++)
		{
			Ft=value(ft(1)(i)); //hardwired now just to use first gear (commercial fishery in Pac hake assessment
			equilibrium(Ft,value(ro),value(kappa),value(m_bar),age,wt_obs(i),fec(i),va_bar,re,ye,be,phiq,dphiq_df,dre_df,dspr_df, dsb_df, spre, ue);
			ispr(i) = (-ro/((kappa-1)*re-ro*kappa));
		}
	}
	//5. Run equilibrium again, similar to first run, to find f that reduces SPR to 0.4
	//The resultin f40 from this have been checked against the ccam.eql table produced at step 3 above - i.e. this f40 gives SPR=0.4
	dvariable F40=value(m_bar); //initialise
	double stp=0.001;
	for(i = 1; i<=20; i++)
	{
		equilibrium(value(F40),value(ro),value(kappa),value(m_bar),age,avg_wt,avg_fec,va_bar,re,spr40_y,spr40_sb,phiq,dphiq_df,dre_df,dspr_df, dsb_df, spr40_spr, spr40_u);
		spr4 = value(-ro/((kappa-1)*re-ro*kappa)); //spawning potential ratio (should be phie/phie0)
		F40 -= (spr4-0.4)/dspr_df;
		if(sfabs(spr4-0.4)<=1.e-6)  break;
	}
	f40=F40;
	//6. Finally, run equilibrium again, similar to first run, to find f that reduces spawning biomass depletion to 0.4
	//The resultin fb40 from this have been checked against the ccam.eql table produced at step 3 above - i.e. this f40 gives SPR=0.4
	dvariable FB40=value(m_bar); //initialise
	double dde;
	for(i = 1; i<=30; i++)
	{
		equilibrium(value(FB40),value(ro),value(kappa),value(m_bar),age,avg_wt,avg_fec,va_bar,re,sb40_y,sb40_sb,phiq,dphiq_df,dre_df,dspr_df, dsb_df, sb40_spr, sb40_u);
		dde = sb40_sb/value(bo); //CG RF change - added this calculation
		FB40 -= (dde-0.4)/dsb_df;
		//cout<<FB40<<"  "<<dde<<endl;
		if(sfabs(dde-0.4)<=1.e-6)  break;
	}
	fb40=FB40;
	sb40_de=dde;

	if(verbose)cout<<"**** Ok after calc_reference_points ****"<<endl;
  }

FUNCTION void simulation_model(const long& seed)
  {
	/*
	Call this routine to simulate data for simulation testing.
	The random number seed can be used to repeat the same
	sequence of random number for simulation testing.

	Implemented using the "-SimFlag 99" command line option where
	99 is the random number seed.

	-SimFlag 99 is a special case used for the manuscript for case 1.
	-SimFlag 000 is a special case with 0 error (exact data)

	-This routine will over-write the observations in memory
	with simulated data, where the true parameter values are
	the initial values.  Change the standard deviations of the
	random number vectors epsilon (observation error) or
	recruitment devs wt (process error).
	*/

	cout<<"___________________________________________________\n"<<endl;
	cout<<"  **Implementing Simulation--Estimation trial**    "<<endl;
	cout<<"___________________________________________________"<<endl;
	cout<<"\tRandom Seed No.:\t"<< rseed<<endl;
	cout<<"___________________________________________________\n"<<endl;
	//Indexes:
	int i,j,k,ii,ki;

	//3darray Chat(1,ngear,syr,nyr,sage,nage);
	//C.initialize();

	/*----------------------------------*/
	/*	-- Generate random numbers --	*/
	/*----------------------------------*/
	random_number_generator rng(seed);
	dvector wt(syr-nage-1,nyr);			//recruitment anomalies
	dmatrix epsilon(1,nit,1,nit_nobs);  //observation errors in survey

	double sig = value(rho/varphi);
	double tau = value((1.-rho)/varphi);

	if(seed==000)
	{
		cout<<"No Error\n";
		sig=0;
		tau=0;
	}
	wt.fill_randn(rng); wt *= tau;
	epsilon.fill_randn(rng);

	//now loop over surveys and scale the observation errors
	for(k=1;k<=nit;k++)
	{
		for(j=1;j<=nit_nobs(k);j++)
			epsilon(k,j) *= sig/it_wt(k,j);
	}

	/*----------------------------------*/
  /*		--Initialize model--	      	*/
	/*CHANGED now calculating phie      */
  /*        based on m_bar and avg_fec*/
	/*----------------------------------*/
	dvector lx=pow(exp(-value(m_bar)),age-min(age));
	lx(nage)/=(1.-exp(-value(m_bar)));
	double phie=(lx*exp(-value(m_bar)*cntrl(13)))*avg_fec;//fec(syr);
	so=kappa/phie;
	if(cntrl(2)==1) beta=(kappa-1.)/(ro*phie);
	if(cntrl(2)==2) beta=log(kappa)/(ro*phie);
	//Initial numbers-at-age with recruitment devs
	/*for(i=syr;i < syr+sage;i++)
			N(i,sage)=exp(log_avgrec+wt(i));
		for(j=sage+1;j<=nage;j++)
			N(syr,j)=exp(log_avgrec+wt(syr-j))*lx(j);
		*/
	N.initialize();
	if(cntrl(5)){	//If initializing in at unfished conditions
		log_rt(syr) = log(ro);
		for(j=sage;j<=nage;j++)
		{
			N(syr,j)=ro*exp(-m_bar*(j-1.));
		}
	}
	else{			//If starting at unfished conditions
		log_rt(syr) = log_avgrec;
		N(syr,sage)=mfexp(log_rt(syr));
		for(j=sage+1;j<=nage;j++)
		{
			N(syr,j)=mfexp(log_recinit+init_log_rec_devs(j))*exp(-m_bar*(j-sage));
		}
	}
	N(syr,nage)/=(1.-exp(-m_bar));

	//log_rt=log_avgrec+log_rec_devs;
	//log_rt(syr) = log(ro);

	for(i=syr+1;i<=nyr;i++){
		log_rt(i)=log_avgrec+log_rec_devs(i);
		N(i,sage)=mfexp(log_rt(i));
	}
	N(nyr+1,sage)=mfexp(log_avgrec);

	/*
	for(j=sage;j<=nage;j++)
	{
		if(cntrl(5))  //if starting at unfished state
		{
			N(syr,j)=ro*exp(-m_bar*(j-1));
		}
		else{
			log_rt(syr-j+sage)=log_avgrec+log_rec_devs(syr-j+sage);
			N(syr,j)=mfexp(log_rt(syr-j+sage))*exp(-m_bar*(j-sage));
		}
	}

	N(syr,nage)/=(1.-exp(-m_bar));
	*/

	/*----------------------------------*/
    /*		--    Selectivity   --		*/
	/*----------------------------------*/
	/*
		-Based on values in the control file.
		-Add autocorrelated random numbers
		for time varying or something to that
		effect.

		-If seed==99 then set up a special case
		for the cubic spline manunscript using
		the eplogistic function where g goes from
		strongly domed to asymptotic, e.g.,
		g = 0.2 * (nyr-i)/(nyr-syr);
	*/

	/*CHANGED May 15, 2011 calcSelectivities gets called from PRELIMINARY_CALCS*/
	dmatrix va(1,ngear,sage,nage);			//fishery selectivity
	d3_array dlog_sel(1,ngear,syr,nyr,sage,nage);
	dlog_sel=value(log_sel);
	/*
	for(k=1;k<=ngear;k++)
		for(i=syr;i<=nyr;i++)
		{
			//sel(k)(i)=plogis(age,ahat(k),ghat(k));
			log_sel(k)(i)=log(plogis(age,ahat(k),ghat(k)));
			log_sel(k)(i) -= log(mean(exp(log_sel(k)(i))));
		}
		//log_sel(j)(i) -= log(mean(mfexp(log_sel(j)(i))));
	*/

	/*----------------------------------*/
    /*	--  Population dynamics  --		*/
	/*----------------------------------*/
	dmatrix zt(syr,nyr,sage,nage);			//total mortality
	zt.initialize();
	dmatrix ft(syr,nyr,1,ngear);
	ft.initialize();
	dvector sbt(syr,nyr+1);
	sbt.initialize();

	for(i=syr;i<=nyr;i++)
	{
		//total biomass at age
		//dvector bt = elem_prod(value(N(i)),wa);
		dvector bt = elem_prod(value(N(i)),wt_obs(i));

		/*calculate instantaneous fishing mortalities
		based on Baranov's catch equation and the
		observed catch from each fleet.*/
		dvector oct = trans(obs_ct)(i);

		for(k=1;k<=ngear;k++)
			va(k)=exp(dlog_sel(k)(i));
		//get_ft is defined in the Baranov.cxx file
		//CHANGED these ft are based on biomass at age, should be numbers at age
		//ft(i) = get_ft(oct,value(m),va,bt);
		ft(i) = get_ft(oct,value(m),va,value(N(i)),wt_obs(i));
		// overwrite observed catch incase it was modified by get_ft
		for(k=1;k<=ngear;k++)
			obs_ct(k,i)=oct(k);
		//total age-specific mortality
		//dvector zt(sage,nage);
		zt(i)=value(m);
		for(k=1;k<=ngear;k++){
			zt(i)+= ft(i,k)*exp(dlog_sel(k)(i));
		}
		//CHANGED definition of spawning biomass based on ctrl(12)
		sbt(i) = value(elem_prod(N(i),exp(-zt(i)*cntrl(13)))*fec(i));
		//Update numbers at age
		if(i>=syr+sage-1)
		{
			double rt;
			//double et=value(N(i-sage+1))*fec(i-sage+1);
			double et=sbt(i-sage+1);
			if(cntrl(2)==1)rt=value(so*et/(1.+beta*et));
			if(cntrl(2)==2)rt=value(so*et*exp(-beta*et));
			N(i+1,sage)=rt*exp(wt(i)-0.5*tau*tau);
			/*CHANGED The recruitment calculation above is incosistent
			  with the assessment model.  Below recruitment is based on
			  rt=exp(log_avgrec + wt + rt_dev), where the rt_dev calculation
			is based on the BH or Ricker model.*/
			//double rt_dev = log(rt)-value(log_avgrec);
			//N(i+1,sage)=exp(log_avgrec+wt(i));
		}
		N(i+1)(sage+1,nage)=++elem_prod(N(i)(sage,nage-1),exp(-zt(i)(sage,nage-1)));
		N(i+1,nage)+=N(i,nage)*exp(-zt(i,nage));
		//Catch & Catch-at-age
		for(k=1;k<=ngear;k++)
		{
			if(ft(i,k)>0)
			{
				dvector sel = exp(dlog_sel(k)(i));
				d3C(k)(i)=elem_prod(elem_div(ft(i,k)*sel,zt(i)),elem_prod(1.-exp(-zt(i)),value(N(i))));
				obs_ct(k,i)=d3C(k)(i)*wt_obs(i);
			}
			else	//if this is a survey
			{
				dvector sel = exp(dlog_sel(k)(i));
				d3C(k)(i)=elem_prod(elem_div(sel,zt(i)),elem_prod(1.-exp(-zt(i)),value(N(i))));
			}
		}
	}
	//initial values of log_ft_pars set to true values
	ki=1;
	for(k=1;k<=ngear;k++)
		for(i=syr;i<=nyr;i++)
			if(obs_ct(k,i)>0){
				log_ft_pars(ki++)=log(ft(i,k));
			}
	// Error handler to inform user population went extinct.
	if(min(sbt(syr,nyr))<=1.e-5)
	{
		cout<<"---------------------------------------\n";
		cout<<"Simulated population went extinct, try\n";
		cout<<"increasing steepness, Ro and Rbar\n";
		cout<<sbt<<endl;
		cout<<"Minimum spawning biomass="<<min(sbt(syr,nyr))<<endl;
		cout<<"---------------------------------------\n";
		exit(1);
	}
	//Average recruitment calculation
	cout<<"	log(mean(column(N,sage))) = "<<mean(log(column(N,sage)))<<endl;
	cout<<"	log_avgrec = "<<log_avgrec<<endl;

	/*----------------------------------*/
    /*	--  Observation models  --		*/
	/*----------------------------------*/
	//Simulated Age-compositions
	int ig;
	for(k=1;k<=na_gears;k++)
	{
		for(i=1;i<=na_nobs(k);i++)
		{
			ii=A(k,i,a_sage(k)-2);	//index for year
			ig=A(k,i,a_sage(k)-1);	//index for gear
			dvector pa = d3C(ig)(ii);	//
			pa/=sum(pa);

			dvector t1=pa(a_sage(k),a_nage(k));
			t1/=sum(t1);
			A(k)(i)(a_sage(k),a_nage(k))=rmvlogistic(t1,0.3,i+seed);
			if(seed==000)
			{
				A(k)(i)(a_sage(k),a_nage(k))=t1;
			}
		}
	}

	//Relative abundance indices
	//CHANGED fixed this to reflect survey timing etc & survey_type
	for(k=1;k<=nit;k++)
	{
		for(i=1;i<=nit_nobs(k);i++)
		{
			ii=iyr(k,i);
			ig=igr(k,i);
			dvector sel = exp(dlog_sel(ig)(ii));
			dvector Np = value(elem_prod(N(ii),exp(-zt(ii)*it_timing(k,i))));
			switch(survey_type(k))
			{
				case 1: //survey based on numbers
					Np = elem_prod(Np,sel);
				break;
				case 2: //survey based on biomass
					Np = elem_prod(elem_prod(Np,sel),wt_obs(ii));
				break;
				case 3: //survey based on spawning biomass
					Np = elem_prod(Np,fec(ii));
				break;
			}
			it(k,i) = sum(Np) * exp(epsilon(k,i));
		}
	}
	//CHANGED Fixed bug in reference points calc call from simulation model,
	//had to calculate m_bar before running this routine.
	calc_reference_points();
	cout<<"___________________________________________________"<<endl;
	ofstream ofs("ccam.sim");
	ofs<<"fmsy\n"<<fmsy<<endl;
	ofs<<"msy\n"<<msy<<endl;
	ofs<<"bmsy\n"<<bmsy<<endl;
	ofs<<"va\n"<<va<<endl;
	ofs<<"sbt\n"<<sbt<<endl;//<<rowsum(elem_prod(N,fec))<<endl;
	ofs<<"rt\n"<<rt<<endl;
	ofs<<"ct\n"<<obs_ct<<endl;
	ofs<<"ft\n"<<trans(ft)<<endl;
	ofs<<"ut\n"<<elem_div(colsum(obs_ct),N.sub(syr,nyr)*wa)<<endl;
	ofs<<"iyr\n"<<iyr<<endl;
	ofs<<"it\n"<<it<<endl;
	ofs<<"N\n"<<N<<endl;
	ofs<<"A\n"<<A<<endl;
	ofs<<"dlog_sel\n"<<dlog_sel<<endl;
	cout<<"  -- Simuation results written to ccam.sim --\n";
	cout<<"___________________________________________________"<<endl;
  }

FUNCTION dvector cis(const dvector& na)
  {
	//Cohort Influenced Selectivity
	//This function returns a vector of residuals from a
	//linear regression of log(pa)= a+b*age+res that can be
	//used to modify age-based selectivity according to relative
	//cohort strengths.
	//SM  Currently not used at all in ccam and should be deprecated.
	dvector y = log(na);
	dvector x = age;
	double b = sum(elem_prod(x-mean(x),y-mean(y)))/sum(square(x-mean(x)));
	double a = mean(y)-b*mean(x);
	dvector res = y - (a+b*x);
	return(res);
  }

REPORT_SECTION
  {
	if(verbose)
		cout<<"Start of Report Section..."<<endl;
	//cout<<"Female sbo report =  "<<bo/2<<endl<<endl;
	//CG RF change - a large number of changes to report section
	int i,j,k,g;

  // Parameterization selection (management vs biological) and Priors
  // output
	report<<"ManagementOriented\n"<<useManagementParameters<<endl;
	report<<"npar\n"<<npar<<endl;
	REPORT(theta_control)

	report<<"nf\n"<<nf<<endl;
	report<<"runtime\n"<<(long(difftime(finish,start))%3600)%60<<endl;
	report<<"Likelihood\n"<<Lf<<endl;
	report<<"priors\n"<<Pf<<endl;
	report<<"ObjectiveFunction\n"<<f<<endl;

  report<<"MaxGrad\n"<<objective_function_value::gmax<<endl;
  report<<"ExitCode\n"<<iexit<<endl;
  report<<"HangCode\n"<<ihang<<endl;

	REPORT(ControlFile);
	REPORT(f);
	REPORT(nlvec);
	REPORT(pvec);
	REPORT(qvec);
	REPORT(priors);
  // report q priors (always on)
	REPORT(q_prior);
	REPORT(q_mu);
	REPORT(q_sd);
	// report ghat priors (can be on or off)
	REPORT(ghat_p1);
	REPORT(ghat_p2);
	REPORT(ghat_pswitch);

	report<<"\n#______________________________\n"<<endl;
	report<<"#______________________________"<<endl;
	report<<"#  Estimated parameters"<<endl;
	if(useManagementParameters) REPORT(msy);
	if(useManagementParameters) REPORT(fmsy);
	REPORT(theta);  // all 7 parameter estimates.
	REPORT(h);
	REPORT(ro);
	REPORT(m);
	REPORT(log_avgrec);
	double rbar=value(exp(log_avgrec));
	REPORT(rbar);
	REPORT(log_recinit);
	double rinit=value(exp(log_recinit));
	REPORT(rinit);
	double tau = value((1.-rho)/varphi);
	double sig = value(rho/varphi);
	REPORT(tau);
	REPORT(sig);
	REPORT(age_tau2);
	REPORT(rho);
	REPORT(varphi);
	REPORT(q);
	REPORT(qt);
	REPORT(so); //SR alpha parameter
	//REPORT(beta);//SR beta parameter
	report<<"bta\n"<<beta<<endl; //female spawning biomass
	report<<"sbo\n"<<bo/2.<<endl; //female spawning biomass
	REPORT(kappa); //recruitment compensation ratio
	//Selectivity
	report<<"#______________________________"<<endl;
	report<<"#  Selectivity parameters"<<endl;

	report<<"ahat"<<endl;
	for(k=1;k<=ngear;k++) report<<mfexp(sel_par(k,1,1))<<" ";
	report<<endl<<"ghat"<<endl;
	for(k=1;k<=ngear;k++) report<<mfexp(sel_par(k,1,2))<<" ";
	report<<endl;

	report<<"log_sel"<<endl;
	for(k=1;k<=ngear;k++)
		for(i=syr;i<=nyr;i++)
			report<<k<<"\t"<<log_sel(k)(i)<<endl;
	report<<"selectivity"<<endl;
	for(g=1;g<=ngear;g++){
		for(i=syr;i<=nyr;i++) report<<g<<"\t"<<mfexp(log_sel(g)(i))<<endl;
	}
	report<<"#______________________________"<<endl;
	report<<"#  Data and other parameters and variables"<<endl;

	ivector yr(syr,nyr);
	ivector yrs(syr,nyr+1);
	yr.fill_seqadd(syr,1); \
	yrs.fill_seqadd(syr,1); 
	REPORT(ngear);
	REPORT(nits);
	REPORT(yr);
	REPORT(yrs);
	REPORT(iyr); //WHEN MULTIPLE SURVEYS - MAKE SURE THAT CODE IS ROBUST TO FIRST SURVEY BEING SHORTER THAN SUBSEQUENT SURVEYS
	REPORT(age);
	REPORT(la);
	REPORT(wa);
	REPORT(fec);
	REPORT(obs_ct);
	REPORT(ct);
	REPORT(ft);
	/*FIXED small problem here with array bounds if using -retro option*/
	report<<"ut\n"<<elem_div(colsum(obs_ct)(syr,nyr),N.sub(syr,nyr)*wa)<<endl;
	int rectype=int(cntrl(2));
	REPORT(rectype);
	REPORT(rt);
	dvector ln_rt=value(log_rt(syr,nyr));
	REPORT(ln_rt);
	REPORT(delta);
	REPORT(it); //WHEN MULTIPLE SURVEYS - MAKE SURE THAT CODE IS ROBUST TO FIRST SURVEY BEING SHORTER THAN SUBSEQUENT SURVEYS
	REPORT(pit); //WHEN MULTIPLE SURVEYS - MAKE SURE THAT CODE IS ROBUST TO FIRST SURVEY BEING SHORTER THAN SUBSEQUENT SURVEYS
	REPORT(epsilon);
	REPORT(F);
	REPORT(M_tot);
	REPORT(a_sage);
	REPORT(a_nage);
	REPORT(A); 
	REPORT(Ahat);
	REPORT(A_nu); //array of residuals
	REPORT(N);
	REPORT(wt_obs);
	//biomasses
	cout<<"Report Section CHECK ***********************************..."<<endl;
	cout<<"N min row index: "<<N.rowmin()<<endl;
	cout<<"N max row index: "<<N.rowmax()<<endl;
	cout<<"N min col index: "<<N.colmin()<<endl;
	cout<<"N max col index: "<<N.colmax()<<endl;
	cout<<"wt_obs min row index: "<<wt_obs.rowmin()<<endl;
	cout<<"wt_obs max row index: "<<wt_obs.rowmax()<<endl;
	cout<<"wt_obs min col index: "<<wt_obs.colmin()<<endl;
	cout<<"wt_obs max col index: "<<wt_obs.colmax()<<endl<<endl;
	dmatrix nt3=value(trans(N).sub(3,nage));
	dmatrix wt3=value(trans(wt_obs).sub(3,nage));
	dvector bt3=colsum(elem_prod(nt3,wt3));

	report<<"bt3\n"<<bt3<<endl;
	report<<"Btotal\n"<<rowsum(elem_prod(N,wt_obs))<<endl;
	dmatrix nt2=value(trans(N).sub(2,nage));
	dmatrix wt2=value(trans(wt_obs).sub(2,nage));
	dvector bt2=colsum(elem_prod(nt2,wt2));
	report<<"bt2\n"<<bt2<<endl;
	dvector rt=value(column(N.sub(syr+1,nyr+1),1)/exp(-m));
	report<<"age0.rt\n"<<rt<<endl;
	dvector rt3(1,3);
	report<<"sbt\n"<<sbt/2.<<endl; //FEMALE SPAWNING BIOMASS

	if(last_phase())
	{
		dvector rt3 = age3_recruitment(value(column(N,3)),wt_obs(nyr+1,3),value(M_tot(nyr,3)));
		REPORT(rt3);
	}

	dvector future_bt = value(elem_prod(elem_prod(N(nyr+1),exp(-M_tot(nyr))),wt_obs(nyr+1)));
	REPORT(future_bt);
	double future_bt4 = sum(future_bt(4,nage));
	REPORT(future_bt4);

	//Parameter controls
	dmatrix ctrl=theta_control;
	REPORT(ctrl);
	REPORT(cntrl);
	//if(last_phase()) projection_model(0);

	if(last_phase())
	{	calc_reference_points();
		if(!useManagementParameters) REPORT(msy);
		if(!useManagementParameters) REPORT(fmsy);
		REPORT(f40); //fishing mortality that reduces spr to 0.4
		REPORT(fb40); //fishing mortality that reduces spr to 0.4
		report<<"bmsy\n"<<bmsy/2.<<endl;
		//REPORT(bmsy);
		report<<"oneminus_ispr\n"<<(1.-ispr)<<endl;
		report<<"oneminus_spr_msy\n"<<(1.-spr_msy)<<endl;
		report<<"sprmsy_status\n"<<(1.-ispr)/(1.-spr_msy)<<endl;
		report<<"sbtmsy_status\n"<<sbt/bmsy<<endl;
		//REPORT(ispr); //spr implied by current f
		//REPORT(spr_msy); //spr at fmsy (This is the reference point for fmsy based reference points). 
		//Note that, by definition, 0.4 is reference point for f40 based reference point, since spr is already a ratio - for the unfished population spr =1.

		//run forecasting model (mle case) - mcmc case is run in mcmc function
		ofstream ofml("ccam_forecast_mle.rep"); //mle estimates
  	ofml<<"Year\t"<<"CtStream\t"<<"Ft\t"<<"ABC\t"<<"OY\t"<<"SBt\t"<<"depletion\t"<<"SPR40status\t"<<"SPRfmsystatus"<<endl;
		dvector ct_stream(nyr+1,pyr);  //pyr set in data section (from control file value)
		dvector oyt_stream(nyr+1,pyr); //this is the optimal yield - updated by reference by forecast model - these will be passed back into the function for two cases
		//Note ct_stream IS used in forecast_model but IS NOT updated
		//oyt_stream IS NOT used in forecast_model but IS updated so catch stream can be set to its values
		ct_stream.initialize();
		oyt_stream.initialize();
		ct_stream=-999;
		forecast_model(ct_stream,oyt_stream); //run under fmsy
		/* //only difference between this and previous is anomalies
		ct_stream=-9991;
		forecast_model(ct_stream,oyt_stream); //run under OY implied by fmsy
		*/
		ct_stream=-888;
		forecast_model(ct_stream,oyt_stream); //run under f40
		/* //only difference between this and previous is anomalies
		ct_stream=-8881;
		forecast_model(ct_stream,oyt_stream); //run under OY implied by f40
		*/
	}
	if(verbose)cout<<"END of Report Section..."<<endl;
  }

FUNCTION mcmc_output
  {
	int i;
	int g;
	adstring tt="\t";
	if(nf==1){
		adstring str_q;
		str_q="lnq";
		ofstream ofs("ccam.mcmc");
		if(useManagementParameters){
			ofs<<"msy fmsy log.m log.rbar log.rinit rho varphi ";
			
			for(g=1;g<=ngear;g++){
			     if(isel_type(g)==1 || isel_type(g)==6 || isel_type(g)==7){
				cout<<"g "<<g<<endl;
				ofs<<"ahat"<<g<<" "<<"ghat"<<g<<" ";
		       }
		   }
			for(i=1;i<=nits;i++){
			  ofs<<"q"<<i<<" ";
		       }
			ofs<<"bo bmsy log.ro h f40 fb40 ";
		}else{
			ofs<<"log.ro h log.m log.rbar log.rinit rho varphi ";
			for(g=1;g<=ngear;g++){
			     if(isel_type(g)==1 || isel_type(g)==6 || isel_type(g)==7){
				ofs<<"ahat"<<g<<" "<<"ghat"<<g<<" ";
		       }
		   }
			for(i=1;i<=nits;i++){
				ofs<<"q"<<i<<" ";
		       }
			ofs<<"bo bmsy msy fmsy f40 fb40 ";
		}
		ofs<<"SSB Age-4 Poor Average Good ";
		for(i=1;i<=nit;i++)ofs<<str_q<<i<<" ";
		ofs<<"f"<<endl;
		ofstream of1("sbt.mcmc");
		ofstream of2("rt.mcmc");
		ofstream of3("1-spr.mcmc");
		ofstream of4("ft.mcmc");
		ofstream of5("bt3.mcmc");
		ofstream of6("bt.mcmc");
		ofstream of7("surveyres.mcmc");
		ofstream of8("sprmsy_status.mcmc");
		ofstream of9("spr40_status.mcmc");
		ofstream of10("sbtdepletion.mcmc");
		ofstream of11("table_i_summary.mcmc");
		ofstream of12("recDevs.mcmc");
		
		//RF NOTE: these don't need to be here but for now I am having trouble with file.exists in R
		//so for now write the files and they can be read into R even if they are empty
		//turn these off once Chris has helped with the path in file.exists in ccamLoadScenarios
		//Help Chris!! :o)
		ofstream of13("surveyres2.mcmc"); //only use if number of surveys (nits) > 1
		ofstream of14("surveyres3.mcmc"); //only use if number of surveys (nits) > 2
		ofstream of15("surveyres4.mcmc"); //only use if number of surveys (nits) == 4 
		of11<<"Sbo"<<tt<<"Ro"<<tt
			<<"SB40.SB"<<tt<<"SB40.SPR"<<tt<<"SB40.U"<<tt<<"SB40.Y"<<tt<<"SB40.Depletion"<<tt
			<<"SPR40.SB"<<tt<<"SPR40.SPR"<<tt<<"SPR40.U"<<tt<<"SPR40.Y"<<tt
			 <<"MSY.SB"<<tt<<"MSY.SPR"<<tt<<"MSY.U"<<tt<<"MSY.Y"<<endl;
	}
	// leading parameters & reference points
	calc_reference_points();
	//cout<<msy<<endl;
	//cout<<"sprmsy"<<endl<<spr_msy<<endl;
	// decision table output
	dvector future_bt = value(elem_prod(elem_prod(N(nyr+1),exp(-M_tot(nyr))),wt_obs(nyr+1)));
	double future_bt4 = sum(future_bt(4,nage));
	double future_bt3 = sum(future_bt(3,nage));
	dvector rt3 = age3_recruitment(value(column(N,3)),wt_obs(nyr+1,3),value(M_tot(nyr,3)));
	dvector bt = value(rowsum(elem_prod(N,wt_obs)));
	dmatrix nt3=value(trans(N).sub(3,nage));
	dmatrix wt3=value(trans(wt_obs).sub(3,nage));
	dvector bt3=colsum(elem_prod(nt3,wt3));
	
	//write mcmc output file
	ofstream ofs("ccam.mcmc",ios::app);
	ofs<<theta;
 	for(g=1;g<=ngear;g++){
	     if(isel_type(g)==1 || isel_type(g)==6 || isel_type(g)==7){
		ofs<<" "<<mfexp(sel_par(g,1,1))<<" "<<mfexp(sel_par(g,1,2))<<" ";
	       }
	   }
 	
	ofs<<" "<<q<<" ";
	
	
	if(useManagementParameters){
		ofs<<" "<<bo/2.<<" "<<bmsy/2.<<" "<<log(ro)<<" "<<h<<" "<<f40<<" "<<fb40<<" ";
	}else{
		ofs<<" "<<bo/2.<<" "<<bmsy/2.<<" "<<msy<<" "<<fmsy<<" "<<f40<<" "<<fb40<<" ";
	}
	ofs<<sbt(nyr)/2.<<" "<<future_bt4<<" "<<future_bt4+rt3<<" ";
	ofs<<log(q)<<" "<<f<<endl;
	// output spawning stock biomass
	ofstream of1("sbt.mcmc",ios::app);
	of1<<sbt(syr,nyr+1)/2.<<endl;
	// output age-1 recruits
	ofstream of2("rt.mcmc",ios::app);
	of2<<rt<<endl;
	// output spr
	ofstream of3("1-spr.mcmc",ios::app);
	of3<<1-ispr<<endl;
	// output ft
	ofstream of4("ft.mcmc",ios::app);
	of4<<ft(1)<<endl;
	// output age-3 biomass
	ofstream of5("bt3.mcmc",ios::app);
	of5<<bt3<<endl;
	// output total biomass
	ofstream of6("bt.mcmc",ios::app);
	of6<<bt<<endl;
	// output survey residuals
	ofstream of7("surveyres.mcmc",ios::app);
	of7<<epsilon(1)<<endl;
	// output spr status (fmsy policy)
	ofstream of8("sprmsy_status.mcmc",ios::app);
	of8<<(1.-ispr)/(1.-spr_msy)<<endl;
	// output spr status (f40 policy)
	ofstream of9("spr40_status.mcmc",ios::app);
	of9<<(1.-ispr)/(1.-0.4)<<endl;
	ofstream of10("sbtdepletion.mcmc",ios::app);
	of10<<sbt(syr,nyr+1)/bo<<endl;
	//Data for summary table in executive summary
	ofstream of11("table_i_summary.mcmc",ios::app);
	of11<<setprecision(4)<<0.5*bo<<tt<<ro<<tt
	<<0.5*sb40_sb<<tt<<sb40_spr<<tt<<sb40_u<<tt<<sb40_y<<tt<<sb40_de<<tt
	<<0.5*spr40_sb<<tt<<spr40_spr<<tt<<spr40_u<<tt<<spr40_y<<tt
	<<0.5*msy_sb<<tt<<msy_spr<<tt<<msy_u<<tt<<msy_y<<endl;
	
	// output recruitment deviations 
	ofstream of12("recDevs.mcmc", ios::app);
	of12<<mfexp(log_rec_devs)<<endl;
	//Write extra files for survey residuals if more than  one survey
	//(can't have more than 4 right now - make dynamic)
	if(nits>1) {
		ofstream of13("surveyres2.mcmc", ios::app);
		of13<<epsilon(2)<<endl;
	}
	if(nits>2) {
		ofstream of14("surveyres3.mcmc", ios::app);
		of14<<epsilon(3)<<endl;
	} 
	if(nits==4) {
		ofstream of15("surveyres4.mcmc", ios::app);
		of15<<epsilon(4)<<endl;
	}
	

  //Run forecast model
	//Need to vectorise the catch stream so different values can be used for different years (in the case of re-running under OY)
	dvector c_stream(nyr+1,pyr);  //pyr set in data section (from control file value)
	dvector oy_stream(nyr+1,pyr); //this is the optimal yield - updated by reference by forecast model - these will be passed back into the function for two cases
	//Note c_stream IS used in forecast_model but IS NOT updated
	//oy_stream IS NOT used in forecast_model but IS updated so catch stream can be set to its values
	c_stream.initialize();
	oy_stream.initialize();

	//DO NOT CHANGE THE ORDER OF 1-4 BELOW. DO NOT INSERT OTHER SCENARIOS IN THE SERIES
	//THE OY UPDATED IN 1 IS USED IN 2. THE OY UPDATED IN 3 IS USED IN 4.

	//Fmsy scenarios
	//1. Apply Fmsy directly
	c_stream=-999;
	forecast_model(c_stream, oy_stream);
	/*
	//2. Apply OY  -the only difference between this and previous is different anomalies 
	//run again with c_stream set to oy_stream - ie catch rather than F implementation
	c_stream=-9991;
	forecast_model(oy_stream, oy_stream);
	*/
	//F40 scenarios
	//3. Apply F40 directly
	// OY stream has been updated from previous call
	//For the -8881 and -9991 scenarios, forecast_model will read the tac from the oy_stream argument
	c_stream=-888;
	forecast_model(c_stream, oy_stream);
	/*
	//4. Apply OY -the only difference between this and previous is different anomalies
	//run again with c_stream set to oy_stream - ie catch rather than F implementation
	//for(i=nyr+1;i<=pyr; i++) c_stream(i)=oy_stream(i);
	c_stream=-8881;
	forecast_model(oy_stream, oy_stream);
	*/
	//FIXED CATCH SCENARIOS
	c_stream=0.0;
	forecast_model(c_stream, oy_stream);
	c_stream=0.05;
	forecast_model(c_stream, oy_stream);
	c_stream=0.1;
	forecast_model(c_stream, oy_stream);
	c_stream=0.150;
	forecast_model(c_stream, oy_stream);
	c_stream=0.200;
	forecast_model(c_stream, oy_stream);
	c_stream=0.250;
	forecast_model(c_stream, oy_stream);
	c_stream=-777; //OY catch stream from SS
	forecast_model(c_stream, oy_stream);
	/*
	// Projection model.
	// use value(msy) instead of msy (because msy is now a parameter)
	for(int i=0;i<=10;i++)
	{
		double tac = double(i)/10. * 1.5*value(msy); //End change
		projection_model(tac);
	}
	*/
	// Deviance Information Criterion
	/*
		DIC = pd + Dbar, where:
		pd is the effective number of parameters,
		Dbar is the expected (average deviance)
		Dbar = mean(D(theta))
		pd = Dbar - D(mean(theta))
		Agorithm:
			-for each sample compute Dtotal += 2*f;
			-compute a running Dbar = Dtotal/nf;

	//cout<<initial_params::nvarcalc()<<endl;
	int nvar = initial_params::nvarcalc();
	dvector y(1,nvar);
	static dvector sum_y(1,nvar);
	initial_params::xinit(y);
	sum_y = sum_y + y;
	cout<<y(1,3)<<endl;
	cout<<sum_y(1,3)<<endl;
	get_monte_carlo_value(nvar,y);
	if(nf==2)exit(1);
	*/
  }

FUNCTION dvector age3_recruitment(const dvector& rt, const double& wt,const double& M)
  {
	/*
	This routine returns the poor average and good age-3 recruits
	that is used in constructing the decision table for the pacific
	herring fisheries.

	-1) sort the rt vector from small to large
	-2) find the 33rd and 66th percentiles
	-3) take the averge of the 0-33, 34-66, 67-100
	-4) multiply by the average weight
	-5) return the age-3 recruitment biomass
	*/
	dvector s_rt = sort(rt);
	dvector rbar(1,3);
	double idx = floor((nyr-syr+1.)/3.);
	int ix1 = syr+int(idx);
	int ix2 = syr+int(2.*idx);
	rbar(1) = mean(s_rt(syr,ix1));
	rbar(2) = mean(s_rt(ix1+1,ix2));
	rbar(3) = mean(s_rt(ix2+1,nyr));
	rbar = rbar*wt*exp(-M);
	//cout<<rbar<<endl;
	return(rbar);
  }

FUNCTION void projection_model(const double& tac);
  {
	/*
	This routine conducts population projections based on
	the estimated values of theta.  Note that all variables
	in this routine are data type variables.

	Arguments:
	tac is the total allowable catch that must be allocated
	to each gear type based on allocation(k)

	theta(1) = log_ro
	theta(2) = h
	theta(3) = log_m
	theta(4) = log_avgrec
	theta(5) = log_recinit
	theta(6) = rho
	theta(7) = vartheta
	*/
	int i,j,k;
	// pyr is now set in control file and data section
	// int pyr = nyr+2;	//projection year.
	// --derive stock recruitment parameters
	dvector lx(sage,nage); lx=1;
	for(i=sage+1;i<=nage;i++) lx(i)=lx(i-1)*exp(-value(m_bar));
	lx(nage)/=(1.-exp(-value(m_bar)));
	double phib = lx*fec(nyr);//(lx*exp(-value(m_bar))) * fec(nyr);//avg_fec;
	double so = value(kappa)/phib;		//max recruits per spawner
	double beta;
	double bo = value(ro)*phib;  				//unfished spawning biomass
	switch(int(cntrl(2)))
	{
		case 1:
			beta = (value(kappa)-1.)/bo;
		break;
		case 2:
			beta = log(value(kappa))/bo;
		break;
	}
	dvector p_sbt(syr,pyr);
	dvector p_ct(1,ngear);
	dmatrix p_ft(nyr+1,pyr);
	dmatrix p_N(syr,pyr+1,sage,nage);
	dmatrix p_Z(syr,pyr,sage,nage);
	p_N.initialize();
	p_N.sub(syr,nyr+1) = value(N.sub(syr,nyr+1));
	p_sbt(syr,nyr)=value(sbt(syr,nyr));
	p_Z.sub(syr,nyr) = value(Z.sub(syr,nyr));
	//selecticity
	/*CHANGED User to specifies allocation among gear types in data file.*/
	dmatrix va_bar(1,ngear,sage,nage);
	for(k=1;k<=ngear;k++)
	{
		//va_bar+=allocation(k)*value(exp(log_sel(k)(nyr)));
		/*cout<<exp(log_sel(j)(nyr))<<endl;*/
		p_ct(k)   = allocation(k)*tac;
		va_bar(k) = exp(value(log_sel(k)(nyr)));
	}
	for(i = nyr+1; i<=pyr; i++)
	{
		//get_ft is defined in the Baranov.cxx file
		//(wt_obs(nyr+1) is the average wt at age in the last 5 years)
		p_ft(i) = get_ft(p_ct,value(m_bar),va_bar,p_N(i),wt_obs(nyr+1));
		//Calculate mortality
		p_Z(i) = value(m_bar);
		for(k=1;k<=ngear;k++)
		{
			p_Z(i)+=p_ft(i,k)*va_bar(k);
		}
		//Spawning biomass
		p_sbt(i) = elem_prod(p_N(i),exp(-p_Z(i)*cntrl(13)))*fec(nyr);//avg_fec;
		//Age-sage recruits
		double tau = value((1.-rho)/varphi);
		double xx = randn(int(tac)+i)*tau;
		if(i>=syr+sage-1)
		{
			double rt;
			double et=p_sbt(i-sage+1);
			if(cntrl(2)==1)rt=(so*et/(1.+beta*et));
			if(cntrl(2)==2)rt=(so*et*exp(-beta*et));
			p_N(i+1,sage)=rt*exp(xx-0.5*tau*tau);
		}
		//Update numbers at age
		p_N(i+1)(sage+1,nage)=++elem_prod(p_N(i)(sage,nage-1),exp(-p_Z(i)(sage,nage-1)));
		p_N(i+1,nage)+=p_N(i,nage)*exp(-p_Z(i,nage));
	}
	if(nf==1)
	{
		ofstream ofs(BaseFileName + ".proj");
		ofs<<"TAC\t PSC\t PSS\t Ut"<<endl;
	}
	ofstream ofs(BaseFileName + ".proj",ios::app);
	ofs<<tac<<"\t"
	<<0.25*bo/p_sbt(pyr)<<"\t"
	<<p_sbt(pyr-1)/p_sbt(pyr)<<"\t"
	<<(tac/(p_N(pyr-1)(3,nage)*wt_obs(nyr+1)(3,nage)))<<endl;
  }
  //based on projection_model
  //calculates ABC and OY using harvest control rule
  //Can run with optimal F (if tac is -999 or -888) or with fixed constant tac
  //RF FixMe. Make sure weight and fecundity line up with SS3
FUNCTION void forecast_model(const dvector& TAC, dvector& OY); //, dvector& oy : updates oy by reference
	/*
	SPECIAL CASES - RUNS MODEL UNDER OPTIMAL F AND UPDATES TAC
		-999 IMPLIES FMSY
		-888 IMPLIES F40
		-777 IMPLIES SS catch stream

		//the results from the two scenarios below are exactly the same with no process error. Makes no difference whether F or TAC is applied
		//the differences come because of different random process errors between the ABC being calculated and being applied
		//talk to Ian about SS
		-9991 IMPLIES OY UNDER FMSY
		-8881 IMPLIES OY UNDER F40
		Differs from TINSS - F40 and Fmsy values have already been estimated (if in management mode) or calculated in the function equilibrium. 
		By the time this function is called, they are known
	*/
	//End change

	/*
	This routine conducts population projections based on
	the estimated values of theta.  Note that all variables
	in this routine are data type variables.

	Arguments:
	tac is the total allowable catch that must be allocated
	to each gear type based on allocation(k)
	*/

	int i,j,k;
	dvector tac(nyr+1,pyr); tac.initialize();
	dvector ABC(nyr+1,pyr); ABC.initialize();
	double spre,ue;

	tac=TAC;
	//SPECIAL CASE - overwrite tac
	//If tac = -9991 or -8881 the model is being re-run under the previous scenario's OY
	//   -9991 = Fmsy OY
	//   -8881 = F40 OY
	//If tac = -777 the model is being run under the Stock Synthesis ABC catch stream
	for(i = nyr+1; i<=pyr; i++){
		//Read TAC from OY argument from previous scenario
		if(tac(i) == -9991 || tac(i) == -8881) tac(i) = OY(i);
		//Read TAC from SSstream (SS catch stream)
		if(tac(i) == -777) tac(i) = SSstream(i);
	}
	//cout<<tac<<endl;
	dmatrix p_ct(nyr+1,pyr,1,ngear);//CG RF added gear dimension
	dmatrix p_ft(nyr+1,pyr,1,ngear); 
	dmatrix p_N(syr,pyr+1,sage,nage);
	dmatrix p_Z(syr,pyr,sage,nage);
	dvector p_sbt(syr,pyr); //Projected sbt
	dvector prt(nyr+1,pyr); //projected recruits
	dvector pspr(nyr+1,pyr); //Projected spr

	double re,ye,be,phiq,dphiq_df,dre_df,fe;
	double dye_df,ddye_df,spr,dspr_df, dsb_df;

	// --derive stock recruitment parameters
	dvector lx(sage,nage); lx=1.;
	for(i=sage+1;i<=nage;i++) lx(i)=lx(i-1)*exp(-value(m_bar));
	lx(nage)/=(1.-exp(-value(m_bar)));

	double phib = lx*avg_fec; //(lx*exp(-value(m_bar))) * fec(nyr);//;  avg_wt
	double so = value(kappa)/phib;		//max recruits per spawner
	double beta;
	double sbo = value(ro)*phib;  		//unfished spawning biomass RF fixme - this is based on ave fecundity- make consistent through model and with SS3	
	switch(int(cntrl(2)))
	{
		case 1:
			beta = (value(kappa)-1.)/sbo;
		break;
		case 2:
			beta = log(value(kappa))/sbo;
		break;
	}

	//Fill matrices and vectors with historical values
	p_ft.initialize();
	p_ct.initialize();
	p_N.initialize();
	p_N.sub(syr,nyr+1) = value(N.sub(syr,nyr+1));
	p_sbt(syr,nyr+1)=value(sbt(syr,nyr+1));
	p_Z.sub(syr,nyr) = value(Z.sub(syr,nyr));

	//selectivity
	/* User to specifies allocation among gear types in data file.*/
	//RF: for hake all catch is allocated to gear 1 (i.e., commercial fishery)
	dmatrix va_bar(1,ngear,sage,nage);
	for(k=1;k<=ngear;k++)
	{
		//va_bar+=allocation(k)*value(exp(log_sel(k)(nyr)));
		/*cout<<exp(log_sel(j)(nyr))<<endl;*/
		//allocate F if an optimal F is to be used
		for(i=nyr+1; i<=pyr;i++){
			//cout<<"f40 = "<<i<<"  "<<value(f40)<<endl;
			if(tac(i)>=0.) p_ct(i)(k)   = allocation(k)*tac(i); //allocate catch if a catch has been passed in
			if(tac(i)==-999) p_ft(i)(k) = allocation(k)*value(fmsy);
			if(tac(i)==-888) p_ft(i)(k) = allocation(k)*value(f40);
		}	//End change
		va_bar(k) = exp(value(log_sel(k)(nyr)));
	}
	for(i = nyr+1; i<=pyr; i++)
	{
		//get_ft is defined in the Baranov.cxx file
		//(wt_obs(nyr+1) is the average wt at age in the last 5 years)
		//if not forecasting under optimal F, just do straight forecast under catch with no harvest control rule
		ABC(i)=tac(i);
		OY(i)=tac(i);
		//cout<<"p_ft(i)"<<p_ft(i)<<endl;
		//if forecasting under optimal F, get ABC (before applying harvest control rule)
		//then adjust F if necessary using harvest control rule
		//then calculate OY under harvest control rule
		//fixme - use this with tac to check that get_ft is working (i.e. remove the if statement and run for all tac scenarios)
		if(tac(i) == -999 || tac(i) == -888) {
				ABC(i)=get_abc(p_ft(i)(1),value(m_bar), va_bar(1), p_N(i), wt_obs(nyr+1)); //BARANOV
				//Now run harvest control rule (set in control file - currently set to 1 = 40-10 rule)
				// if sbt/sbo > 0.4, p_ft is adjusted down
				p_ft(i)(1)=harvest_control_rule(p_ft(i)(1), hcr, p_sbt(i-1), sbo); //check indexing
				//p_ft(i)(1)=p_ft(i)(1); //test without hcr
				if(verbose){
					cout<<"p_ft(i)2"<<p_ft(i)<<endl;
					cout<<"Check that ft is same as optimal F!! (only when sbt>0.4sbo)"<<endl;
				}
				//Now can get OY
				OY(i)=get_abc(p_ft(i)(1),value(m_bar), va_bar(1), p_N(i), wt_obs(nyr+1));
			}else
			{ p_ft(i) = get_ft(p_ct(i),value(m_bar),va_bar,p_N(i),wt_obs(nyr+1)); }	// NOTE this returns a matrix of F (1 column for each gear - one column will be zero if all catch is allocated to one fleet)

		//Calculate mortality
		p_Z(i) = value(m_bar);
		for(k=1;k<=ngear;k++)
		{
			p_Z(i)+=p_ft(i,k)*va_bar(k);
		}//end for
		//Spawning biomass
		p_sbt(i) = elem_prod(p_N(i),exp(-p_Z(i)*cntrl(13)))*fec(nyr+1);//avg_fec;

		//Age-sage recruits
		double tau = value((1.-rho)/varphi);
		double xx = randn(int(tac(i))+i)*tau;
		if(i>=syr+sage-1)
		{
			double rt;
			double et=p_sbt(i-sage+1);
			if(cntrl(2)==1)rt=(so*et/(1.+beta*et));
			if(cntrl(2)==2)rt=(so*et*exp(-beta*et));
			p_N(i+1,sage)=rt*exp(xx-0.5*tau*tau);
		}
		//Update numbers at age
		p_N(i+1)(sage+1,nage)=++elem_prod(p_N(i)(sage,nage-1),exp(-p_Z(i)(sage,nage-1)));
		p_N(i+1,nage)+=p_N(i,nage)*exp(-p_Z(i,nage));
		equilibrium(p_ft(i)(1),value(ro),value(kappa),value(m_bar),age,wt_obs(nyr+1), fec(nyr+1),va_bar(1),re,ye,be,phiq,dphiq_df,dre_df,dspr_df, dsb_df, spre, ue);
		pspr(i) = (-value(ro)/((value(kappa)-1)*re-value(ro)*value(kappa)));
		//cout<<"sum(p_N(i))"<<sum(p_N(i))<<endl;
		//cout<<"sum(p_N(i+1))"<<sum(p_N(i+1))<<endl;
	}
	//Get Rt_2009 and forecast recruitment estimate for output file (Phake assessment)
	double rt09=value(rt(2009));
	prt=column(p_N.sub(nyr+1,pyr),sage);

	if(nf==1){
		ofstream ofmc("ccam_forecast_mcmc.rep");
		ofmc<<"Year\t"<<"CtStream\t"<<"Ft\t"<<"ABC\t"<<"OY\t"<<"Rt_2009\t"<<"ForecastRt\t"<<"Sbt\t"<<"depletion\t"<<"SPR40status\t"<<"SPRmsystatus\t"<<"SBmsy"<<endl;
	}
	ofstream ofmc("ccam_forecast_mcmc.rep",ios::app);  //appendposterior forecast results (mceval phase)
	ofstream ofml("ccam_forecast_mle.rep",ios::app); //appendmle forecase results
	for(int i=nyr+1;i<=pyr;i++){
		if(!mceval_phase()){
			cout<<i<<"\t"<<TAC(i)<<"\t"<<p_ft(i)(1)<<"\t"<<ABC(i)<<"\t"<<OY(i)<<"\t"<<p_sbt(i)/2.<<"\t"<<p_sbt(i)/sbo<<"\t"<<(1-pspr(i))/(1-0.4)<<"\t"<<(1-pspr(i))/(1-spr_msy)<<endl;
		}

		if(!mceval_phase()){
			ofml<<i<<"\t"<<TAC(i)<<"\t"<<p_ft(i)(1)<<"\t"<<ABC(i)<<"\t"<<OY(i)<<"\t"<<p_sbt(i)/2.<<"\t"<<p_sbt(i)/sbo<<"\t"<<(1-pspr(i))/(1-0.4)<<"\t"<<(1-pspr(i))/(1-spr_msy)<<endl;
		}

		if(mceval_phase()){
			//cout<<"sbo eval = "<<sbo<<endl;
			ofmc<<setprecision(4)<<i<<"\t"<<TAC(i)<<"\t"<<p_ft(i)(1)<<"\t"<<ABC(i)<<"\t"<<OY(i)<<"\t"<<rt09<<"\t"<<prt(i)<<"\t"<<p_sbt(i)/2.<<"\t"<<p_sbt(i)/sbo<<"\t"<<(1-pspr(i))/(1-0.4)<<"\t"<<(1-pspr(i))/(1-spr_msy)<<"\t"<<bmsy/2.<<endl;
		}
	}

FUNCTION double harvest_control_rule(const double& fopt, const double& hcr, const double& sbt, const double& sbo);
  {
  	double foy;
  	if(hcr==1){
		//40:10 Rule
		double target=0.4;
		double limit=0.1;
		double d=sbt/sbo;

		//Calculate 40:10 adjustment to F
		if(d>=target) foy=fopt;
		if(d<=target) foy=4./3.*((d-0.1)/d)*fopt;
		if(d<=limit) foy=0.;
	}
	return(foy);
  }

FUNCTION double get_abc(const double& fopt, const double& m,const dvector& va, const dvector& N, const dvector wa);
	  {	// This is the 40/10 adjustment to the F
		int a, A;
		double abc; //function returns this
		a=N.indexmin(); A=N.indexmax();
		dvector fage(a,A);
		dvector zage(a,A);
		dvector ominus(a,A);
		dvector abc_a(a,A);
		dvector oy_a(a,A);
		dvector bage(a,A);
		bage=elem_prod(N,wa);
		fage=fopt*va;
		zage=fage+m;
		dvector fz=elem_div(fage,zage);
		ominus=(1.-mfexp(-zage));
		abc_a=elem_prod(elem_prod(fz,ominus),bage);
		//cout<<"fopt = "<<fopt<<endl;
		//cout<<"va = "<<va<<endl;
		//cout<<"abc_a = "<<abc_a<<endl;

		/*
		//test line above
		for(int j = a; j<=A; j++)
		{
			abc_a(j)=(fage(j)/zage(j)) * (1.-mfexp(-zage(j)))*bage(j);
			cout<<"abc_a(j) = "<<abc_a(j)<<endl;
			cout<<"bage(j) = "<<bage(j)<<endl;
		}
			double Bt = sum(value(bage));
			cout<<"Bt = "<<Bt<<endl;
			cout<<"sumbage = "<<sum(bage)<<endl;
		 */

		abc=sum(abc_a);
		return(abc);
	}

TOP_OF_MAIN_SECTION
	time(&start);
	arrmblsize = 50000000;
	gradient_structure::set_GRADSTACK_BUFFER_SIZE(4.3e7);//4.5e7
	gradient_structure::set_CMPDIF_BUFFER_SIZE(1.e7);
	gradient_structure::set_MAX_NVAR_OFFSET(5000);
	gradient_structure::set_NUM_DEPENDENT_VARIABLES(5000);

GLOBALS_SECTION
	/**
	\def REPORT(object)
	Prints name and value of \a object on ADMB report %ofstream file.
	*/
	#undef REPORT
	#define REPORT(object) report << #object "\n" << object << endl;

	#undef COUT
	#define COUT(object) cout << #object "\n" << object <<endl;

	#include <admodel.h>
	#include <time.h>
	#include <string.h>
	#include "baranov.cxx"
        #include <statsLib.h>//CG version
       //  #include <stats.cxx>//TEMP UNTIL I PUT NEW BUILD ON MY LAPTOP (RF)
	#include "managementoriented.cxx"
	time_t start,finish;
	long hour,minute,second;
	double elapsed_time;
	bool mcmcPhase = 0;
	bool mcmcEvalPhase = 0;
	adstring BaseFileName;
  adstring ReportFileName;
	adstring stripExtension(adstring fileName)
	{
		/*
		This function strips the file extension
		from the fileName argument and returns
		the file name without the extension.
		*/
		const int length = fileName.size();
		for (int i=length; i>=0; --i)
		{
			if (fileName(i)=='.')
			{
				return fileName(1,i-1);
			}
		}
		return fileName;
	}

FINAL_SECTION
	time(&finish);
	elapsed_time=difftime(finish,start);
	hour=long(elapsed_time)/3600;
	minute=long(elapsed_time)%3600/60;
	second=(long(elapsed_time)%3600)%60;
	cout<<endl<<endl<<"*******************************************"<<endl;
	cout<<"--Start time: "<<ctime(&start)<<endl;
	cout<<"--Finish time: "<<ctime(&finish)<<endl;
	cout<<"--Runtime: ";
	cout<<hour<<" hours, "<<minute<<" minutes, "<<second<<" seconds"<<endl;
	cout<<"--Number of function evaluations: "<<nf<<endl;
	cout<<"--Results are saved with the base name:\n"<<"\t"<<BaseFileName<<endl;
	cout<<"*******************************************"<<endl;
